{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
    "\n",
    "np.random.seed(9876789)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS estimation\n",
    "\n",
    "Artificial data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 100\n",
    "x = np.linspace(0, 10, 100)\n",
    "X = np.column_stack((x, x**2))\n",
    "beta = np.array([1, 0.1, 10])\n",
    "e = np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model needs an intercept so we add a column of 1s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = sm.add_constant(X)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       1.000\n",
      "Model:                            OLS   Adj. R-squared:                  1.000\n",
      "Method:                 Least Squares   F-statistic:                 4.020e+06\n",
      "Date:                Tue, 24 Dec 2019   Prob (F-statistic):          2.83e-239\n",
      "Time:                        14:58:59   Log-Likelihood:                -146.51\n",
      "No. Observations:                 100   AIC:                             299.0\n",
      "Df Residuals:                      97   BIC:                             306.8\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.3423      0.313      4.292      0.000       0.722       1.963\n",
      "x1            -0.0402      0.145     -0.278      0.781      -0.327       0.247\n",
      "x2            10.0103      0.014    715.745      0.000       9.982      10.038\n",
      "==============================================================================\n",
      "Omnibus:                        2.042   Durbin-Watson:                   2.274\n",
      "Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.875\n",
      "Skew:                           0.234   Prob(JB):                        0.392\n",
      "Kurtosis:                       2.519   Cond. No.                         144.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "model = sm.OLS(y, X)\n",
    "results = model.fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quantities of interest can be extracted directly from the fitted model. Type ``dir(results)`` for a full list. Here are some examples:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [ 1.34233516 -0.04024948 10.01025357]\n",
      "R2:  0.9999879365025871\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', results.params)\n",
    "print('R2: ', results.rsquared)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS non-linear curve but linear in parameters\n",
    "\n",
    "We simulate artificial data with a non-linear relationship between x and y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "sig = 0.5\n",
    "x = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x, np.sin(x), (x-5)**2, np.ones(nsample)))\n",
    "beta = [0.5, 0.5, -0.02, 5.]\n",
    "\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + sig * np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.933\n",
      "Model:                            OLS   Adj. R-squared:                  0.928\n",
      "Method:                 Least Squares   F-statistic:                     211.8\n",
      "Date:                Tue, 24 Dec 2019   Prob (F-statistic):           6.30e-27\n",
      "Time:                        14:58:59   Log-Likelihood:                -34.438\n",
      "No. Observations:                  50   AIC:                             76.88\n",
      "Df Residuals:                      46   BIC:                             84.52\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.4687      0.026     17.751      0.000       0.416       0.522\n",
      "x2             0.4836      0.104      4.659      0.000       0.275       0.693\n",
      "x3            -0.0174      0.002     -7.507      0.000      -0.022      -0.013\n",
      "const          5.2058      0.171     30.405      0.000       4.861       5.550\n",
      "==============================================================================\n",
      "Omnibus:                        0.655   Durbin-Watson:                   2.896\n",
      "Prob(Omnibus):                  0.721   Jarque-Bera (JB):                0.360\n",
      "Skew:                           0.207   Prob(JB):                        0.835\n",
      "Kurtosis:                       3.026   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "res = sm.OLS(y, X).fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract other quantities of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [ 0.46872448  0.48360119 -0.01740479  5.20584496]\n",
      "Standard errors:  [0.02640602 0.10380518 0.00231847 0.17121765]\n",
      "Predicted values:  [ 4.77072516  5.22213464  5.63620761  5.98658823  6.25643234  6.44117491\n",
      "  6.54928009  6.60085051  6.62432454  6.6518039   6.71377946  6.83412169\n",
      "  7.02615877  7.29048685  7.61487206  7.97626054  8.34456611  8.68761335\n",
      "  8.97642389  9.18997755  9.31866582  9.36587056  9.34740836  9.28893189\n",
      "  9.22171529  9.17751587  9.1833565   9.25708583  9.40444579  9.61812821\n",
      "  9.87897556 10.15912843 10.42660281 10.65054491 10.8063004  10.87946503\n",
      " 10.86825119 10.78378163 10.64826203 10.49133265 10.34519853 10.23933827\n",
      " 10.19566084 10.22490593 10.32487947 10.48081414 10.66779556 10.85485568\n",
      " 11.01006072 11.10575781]\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', res.params)\n",
    "print('Standard errors: ', res.bse)\n",
    "print('Predicted values: ', res.predict())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare the true relationship to OLS predictions. Confidence intervals around the predictions are built using the ``wls_prediction_std`` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xdc1eUXwPHPlwsITtwDt5aZWppakgs1U0vNsLRyVGpa2lDTwpb9stKyNEfLbFjmzJGr0FCcqKGYe2apaG7Ewbz3+/vjiLi9wJ1w3q8XLwUu3+8D6D33eZ7znGOYpolSSimlXM/H3QNQSimlcisNwkoppZSbaBBWSiml3ESDsFJKKeUmGoSVUkopN9EgrJRSSrmJBmGllFLKTTQIK6WUUm6iQVgppZRyEw3CSimllJv4OvsGxYoVMytWrOjs2yillFIeYcOGDSdM0yxuz2OdHoQrVqxITEyMs2+jlFJKeQTDMP6197G6HK2UUkq5iQZhpZRSyk00CCullFJu4vQ94etJTU3l0KFDJCUlueP2ThUQEEDZsmXx8/Nz91CUUkp5OLcE4UOHDlGgQAEqVqyIYRjuGIJTmKbJyZMnOXToEJUqVXL3cJRSSnk4tyxHJyUlUbRo0RwVgAEMw6Bo0aI5coavlFLK8dy2J5zTAnC6nPp9KaWUcjyvSMyaGxtHwxFLqRS+kIYjljI3Ns5h13733XeJioq67uc2bdrEpk2bHHYvpZRS6nJu2RPOjLmxcQyZvYXEVCsAcfGJDJm9BYAOdYKdeu/0AFy7dm2n3kcppVTu5PFBeGTErksBOF1iqpWREbuyHIRPnz7N448/jtVqxTRN6tWrR+vWrTl//jxVq1bl+++/Z8iQIcyZMweAn376icjISM6dO8djjz12xeOUUkqprPL45ejD8YmZ+rg9JkyYQNu2bVm2bBl+fn4cOXKEl156iT/++IN//vmHo0ePMnz4cMLDwwkPDycyMhLguo9TSimlssrjg3CZoMBMfdwe+/fv5+677wagXr16+Pn5MXHiRLp06cKpU6dITLx+gLf3cUoppbzEmTNuvb3HB+HBraoR6Ge54mOBfhYGt6qW5WuWL1+ebdu2AbLv++233/LYY48xdepU8uXLl3GfwEAuXLgAyBngGz1OKaWUl/r8c3BjkyGPD8Id6gQzPKwWwUGBGEBwUCDDw2plKymrd+/ezJo1i9DQUBISEmjZsiXDhw+nefPmAMTFSfZ1y5YtmT17Ng0bNmTlypU3fJxSSikvsWsX9OgBs2bJ+/36QdGibhuOYZqmU29Qr1498+pWhjt27KB69epOva875fTvTymlvE5sLAwfDr/8AnnywAcfwMCBTrmVYRgbTNOsZ89jPT47WimllMqWfv3giy+gYEEID4f+/aFECXePCrAzCBuGURL4xTTNxoZhlAd+BGzAXqCP6ezptFJKKZUZqang4wMWCzRtCmXLQt++UKiQu0d2hVvuCRuGURiYBKRnIvUBXjBNszlQDqjlvOEppZRSmbRvHzRqBCNGyPudOsGQIR4XgMG+xCwr0BlIADBN803TNHdc/FxR4ISTxqaUUkrZzzThp5+gdm1JwLr9dneP6JZuGYRN00wwTfOag1SGYXQGtpmmedgpI1NKKaXsdeYMdOkC3btDnTqweTM8/ri7R3VLWTqiZBhGZWAQ0P8Gn+9tGEaMYRgxx48fz874nGb8+PGEhoYSGBhIaGjopRKVSimlvNCuXTBnDgwbBsuWQfny7h6RXew+omQYRpRpmqEX94h/B3qZprnlVl/n6UeUqlatyt69ex16TU/6/pRSGebGxjEyYheH4xMpExTI4FbVMldzYN8++PVX2LYN2raFRx+FtDSYOxeKFJHzpkWKQLFiEJj1qn7KTlYrLFkCrVvL+//9B6VKuXdMOP+IUjhQHhh3sXfuUNM0l2fhOoBkiju6W2Dt2vDZZ5n/utDQUOrXr8/mzZuJiIjg3XffJTQ0lNDQUH744QcAOnXqRPfu3Tl27Bi1atXi888/d+zglVJOkeWObDYbvPOOBNqLlfYoUQKqVJEgfPLktcuePj7w6afyBKecIyFBlp8XLICNG2UJ2gMCcGbZvRxtmmboxT9fN02ztGmaoRffshyAPc3atWsJCQkhIiLiho+ZMGECNWvWZMWKFRw5coTNmze7cIRKqay6WUe2KyQnw2+/STlDkIAaGSmBd/RomQ0fPQpvvCGfL1JE9h+XLZMqTBMmwKBBkp0LMssYMkS+TjnG33/D/fdn/J7q1HH3iLLM7cU6sjJjdZaaNWsSFhZ23c8lJiYSGBjIrl27WLNmDVFRUcTHxxMXF8ddd93l4pEqpTLLro5sy5bBs8/Cv//KknLv3uDnBytXgu8Nni79/KDWTU5qrl4NH38sx2WaNYPnnpMZdEBANr6bXCwqCh57TFYoIiKgRQt3jyhbPL52tCvlz5//ivf9/f1JTyz7/fffAahWrRr9+/cnKiqK999/n/JesvmvVG53045sFy7I0nHz5uDvL/u+Bw9KgIUbB2B79OsHBw7A++/D/v3w1FNQo4YEEZV5e/dC8eKwbp3XB2DQIHxT7du3Z9y4cTz//PMUvVjg+7nnnuO3336jSZMmfPXVV5QrV87No1RK2eOmHdmOHoXvvoOXX5bl4/btHTtTDQ6GN9+UJeklS+TvPj5yrnX1asfdJ6dKTYUNG+TvvXpJHejbbnPvmBxEGzg4QU7//pTyVpdnR5fP78socyd133oZDEMCccmSrh3Q/PkS8Nu2lb25KlVce39vcOqUVLyKjoY9e6BMGXeP6JYykx2tM2GlVK7RoU4wq8Obs//Jsiyf+Rp13+mfMRN1dQAGaNUKRo6Ufc4aNSQL+2IPcwVs3w733Sd78l984RUBOLM0CCulcg/TlASp+vXh2DGZiaZnMbuDv79kUu/cCWFhUmiiZUsZZ243e7YE4IQESZh7+ml3j8gp3J4drZRSLjNoEIwaJed6v/zSrc3crxAcDFOmQJ8+MhM2DNkHPXpUuv/kRsuXw513yrGvHPwz0CCslMo9HnoI8uaF996TQOdpmjbN+PuoUdJ4fuRIOdbkkwsWLk+fhiNHJPh+8olUxMrhR7lywW9VKZWrWa1SbAPkSMuwYZ4ZgK/WqRPcey88/zw88IAUqMjJtmyRbYJHHpFSoH5+OT4AQy6dCZ87d45u3bpx/PhxqlSpQrly5bjjjjvo2rXrpcecP3+erl27curUKcqXL8+PP/6I4Q3/cZVSGdLSoFs3mDZNjrXUru3uEdmvUiU5zjRxIrz6qhQE+eYbOWec00yfDj16SL/fWbOydy77MtmuFe4C3jMTjo6G4cPlz2waN24ct912G6tWrSI5OZkZM2Zc85iffvqJkJAQli9fTp48ebj6mJVSysOlpsITT0gAHjHCuwJwOsOQpeht26SQSE47wpSWBoMHy++pTh05CxwS4pBLp9cKj4tPxCSjVvjc2DiHXN9RPGMmHBp67cc6dYK+fSVJoWFDqc1qs8m+yF13wSuvwDPPwIkTUsLsclFRN73dunXr6NWrFwCNGjWidOnS1zwmODiYSZMm8eijjzJx4sSsfV9KKfdITobOnaXy1ahRMGCAu0eUPeXKSSZ3uoEDpVnBgAEZVb28VWysPNePHi3Z4g5ys1rhnjQb9o6Z8JkzGSXebDZ5PxvOnj1Lvnz5AMibNy8JCQnXPKZdu3YMGDCAsLAwXn75ZaxW6zWPUUp5qAULJACPH+/9AfhqVqvUtn79dbj7bli82N0jypzUVEk2O3RIlp0XLpQmDA4MwGBnrXAP4Bkz4ZvNXPPmhZ9/loSKlBT5Rf38c8aSRbFit5z5Xq1gwYKcO3cOkL3fggULXvOYPXv20Lp1azp27EjXrl2ZPHkyT+fQc2pK5TgdO3rfHrC9LBb45RcJXv37S8GP9u1h3DjPb2QfGws9e8qfpgmvvQZ58jjlVmWCAom7TsC9UQ1xd/GOmXBIiGQ3Dhsmf2Zzz+C+++4j6mLgXrlyJYsWLbrmMRMnTmTOnDlYLBZq1qxJUlJStu6plHIyq1XqCqfnjeTEAJzOMKTU5bZtst+9dq27R3RziYnSzrF+fTmCNGuWBGAnummtcA+SK2tHp2dHHz16lNtuu41y5coxZcoUihQpAsAzzzxDWFgYXbp0wTRNChUqxNSpU8mbN69d13f396dUrmOakicybpy8vfiiu0fkWklJcpzHNKWyVMuW0vDeU84WDxkiLxZ69JDzv4ULu+S27sqOzkzt6FwZhJ0tp39/Snmc0aMlWWnAAEnEyq3i4+HBB+HPP2XF8M03oXVrWcJ2tZ074fx5qFtXinBs2CDnnXMBbeCglMo9Zs2Sc7QdO8osKzcLCpKl6e+/h3/+kSXrSpXgr79cc//UVPl9tGgB1atnLDkXLpxrAnBmaRBWSnm3X36BBg3gp588Z/nVnXx85Pjmv//Kz6Zu3Yzeu/PnS0KXM057TJgAFSvKkdF9+6Suw9Spjr9PDuO27GjTNHNkBSpnL+8rpa4yeTKcOweBnpX16nZ+frI60LFjxsc+/VQaI5QrJ0ls7dpJgM6fP/PXP3MGYmLkBVC+fHJ65a674KuvpEa3O5bAvZBbXjYGBARw8uTJHBewTNPk5MmTBOSCeqdKudXx4xJc4uLkyb5QIXePyDssWSKz4+rVYehQuOce6NdPPmea8veRI+WM9fbtcPgwnDwpn9+/X0pm3nefHA0NCpIl5mnT5PP9+sFvv0lg1wBsN7ckZqWmpnLo0KEceewnICCAsmXL4uftVWyU8lSJiVLCcdMm6TPboIG7R+Sd/v1XZrKlSklVwtOnoVo1eYFzuY8+kr3dAweky1PVqlI+s0oVuP12SQTz0lUIZ2VPZyYxyy3L0X5+flSqVMkdt1ZKeTObTRoyrFsHM2dqAM6OChXkLV3hwnDsmGRY79kDu3fLknOTJvL58uVlNpxDpNeWTi9tmV5bGnBpWUvPqJillFL2eO01yb4dNerKvU7lOEFBUlSjfn13j8SpPKW2tAZhpZR3OHcOIiKkEEf//m4Zgje0xlP28ZTa0hqElVLeIX9+WLNG9h/dcLLCU5YvlWOk15a+J24HDQ5sYW35WmwMru7y2tJ6qE4p5dm2bZNyh4mJUKCAwxq+Z9bNli+V9xncqhoh/+1ixs+vM3DlZH6e9iYhR3e7vLa0zoSVUp7r6FF4+GE5g3ryJJQt67aheMrypcqmw4dh0SI69OrFbQWO42PaZDZqTeOt/Meo4eJVDQ3CSinPdOGCtOg7fhxWrHBrAAbvaY2nrhIdLUfZAgNh5UqYN08qhoWGUuOp9vDtGEhJwdffX953MQ3CSinPY7NJN6A//4Q5c6T0opsNblXtij1h8MzWeOoy0dHQrBkkJ8v7QUFSZ7xXLznvXLWqtMeNioLQ0Gy3yc0KDcJKKc+zfz8sXSoNGR55xN2jATKSrzQ72kukpkpwTUuT9318pNPW229f+biQELcE33QahJVSnqdKFSmbWKKEu0dyhQ51gjXoerrkZHj/fZg7V3pL+/tLToG//xWdnExTmktNnw5790rtF3fQIKyU8hxLl8oxpDffhJIl3T2a6zp7Vio+Fismb25K1lbXs3Yt9OwpL+C6d5fa2FctN2/bJoF3+nQpCmaxQMuWkJQE7ij7r/98lFKeYdMmCAuD4GApxpGVzj5OYJqweTP8/rvUCkldsYYm1mUspTnrjBCKFpUJe8mS8meJElA/LZqH8kZRtGOoW5c6c42kJHjjDfjsM0ngW7QI2rSRz4WEsLtoiATe5+TEm4+PxORXX5V/csWKuW/oGoSVUu63e7c0AihYUJ5A3RyAT52ShkPpgffoESv3s4Yhhb6kjVW6BqX5fcDc9t/zwO+vcuFAfs4eyM8Za37MlDRqW2PwJY2k0f6sfT+ShoNC0J4uThIdLSsoERHwwgswYoScJwc2bJDXc6tWyUMbN4bx46XiaalSbhzzZTQIK6Xc68CBjL26JUuubCrgYvv3w+fdovFfE8UyM5S/g+oyPag3DfIvIODcSTjrA0jnOT9bCo9X2QidW1H43Dkpq3nuHIk79uJ3PAULJtiSCXrjBXp88iFVX2hJjz5+lCvntm8vZ7Fa4ZVX4NtvJQnL3x8mToQCBTh1Ct56S1oblygh+X2dO7v9lNt1aRBWSrlXTIycCV6yRFrpuYHVCmPHwrzwNfye0gw/UjHzBGDOi8T37X+gYWvJ0i5cWM4upyf6dOhwxXLz3Ng4po+dwXeTw/GzpmEzfKjis4efTj3M8Q+KMeODzhxo/BSPdoB7k5bj0yz0muXq3FKfOlvf56lT0ts4IkJKmJompKRgWxbF99tDCA+Xh7z8Mvzvf57dbtot/YSVUgrTzKgBfeaM254pN2+WY6Nxf8axNt8DlDu/Uz5hscCwYTBkyJVfEB19w3OlDUcsvaYe8daSVXn02BbePrWVPBHzSEzzw5c0/EnBJ48/PssiL13n6vrUIGeRh4fVcmggdnegz9b3uXGjrCcfPgwDBsirp5QUbL7+9K4SybfbQ2jYED7/HO6+28nfyA1kpp+w1o5WSrne+fOyBP3rr/K+GwJwUpIsWda9x6ThjonsD7yTsql/g5+fBGB/fwm0VwsJkcB8nYSr9BKWG4Or80VIJzYGVyfF148ZZe4h/4Lp+J08Sr7eXQgwUvDFCsmJbOj9FYkXZDLkivrU6QEwLj4Rk4xGFHNj4xx2j1vJ8vf522/QsKGc/V25EkaMIGFOJL/WG0aj5Ejmnwhh0iT5lLsCcGbpcrRSyrWSkyUlNSoKnn/e7i9z5Oxt5Up47jnYtQsGPvoPnyx6ESMkhCX9hzHz91iqbo9h7531eCigPB0ycd1blrYsWBDLM93hp0mYKSkYVht1t/7IqmJHsIwZ5bD61Df7WXlCH90sf5/33COz4NGjoXhxFi2CZ54J4eTJEPq9BIvek6JY3kRnwkop10lLgy5dYPFi+OYbePxxu77MUbO3lBTo1w9Cm1gJObWQ33+HT2dXwli7lrmf/MjLMedYXKgyX4R0YnGhypm+x+BW1Qj0s1zxsWtKW4aEQGQkxrBhGCuWs6ffZ9RMjuHe3ncz9NsZ5IlPu+a6malPfauflSc0orjR93PNx6OjITxcsqrS0uQc2OTJpAYVZ/Bg6e1RurRkQY8d630BGDQIK6VcxTShTx+YNQtGjZL2hHZyxDLthQuSR3Xyi6kcL1CZ74+3pVXelfLJ2rUZuWRPtu/RoU4ww8NqERwUiAEEBwVef58zfUm7cWNuG/8KAQf28Gf9ftQ+sZl/JrXkjpUneCF6BvfE7ch0fepb/azsDoBOZNeLlfS6zx99BDNmwM8/A/DPP3LU6JNP5ETS2rVQu7bLhu5wuhytlHINw5Dzv++8Iwk1mZDd2VtCArRtC/VXfsonDMI4i+z5WjICgaNmiFkpbRkQXJQG68eybVMqYY+v4ds1LxBIIqkWP6In/kLTTFzvVt+HJzSiuGUdbtOUpLj0xgsWCxw+zOzZUhDLZpO4bOdCikfTIKyUcq4TJ6QdYfXqspeXnhGdCdlpI3jiBLRuDbfFzmCk8RpG+oEQqxWWL4f778/2PRylRm0/fn52DbyVjGGCnzWVuz4eAY+3gHz57LrGrb4PT2lEcdMXKwMHShLWxRdJpr8/ozaEMugNqF8fpk2DypVdOFgn0uVopZTz7NkjS6+PPCJ7ej4+WQrCdi1fXkdcHDRtKqUKB71uwadObekre53s56zew9F8moXiE5AH02LBii+Fd0RzMrgWKRHL7Pp6e76PDnWCWR3enP0jHmZ1eHPPO4fcrRuMHAkrVnDilWE8Wy6SQbNCGDhQql/llAAMek5YKeUsq1ZJ8PXxkaNIF2ecWZXZ7Oi//5ZTUIWP7uTThXdIvDVN2US8wTlfd5+fveTiWeTkkFC+HJvKw3N6crTgbZT4cxG3337rL/eY7yMzZs+G9eul7ORFU6ZIGoG/P/zwA7Rr577hZUZmzglrEFZKOd7UqfDMM1CxotSCrlLFpbfftg1aPmAyMP4dXk0djrFmDdx7r0vH4EgLZlzg1T7niEstwaT3/iXMmI2RlOS2RvQOlZYmB7Y/+gjuuw+WLeOCGcgjXc7zx9x85Ak+Rc2uO3irc0XPfyFxUWaCsO4JK+VmXjlruRnThO++kyfUOXOgaFGX3j4mBtq0svFR0iv0SBov5bDq1nXpGBytbae81L4/L127gs+rr2DwKyYGRkAeaV7grYF40SJ46SVZtujTB8aMYfu+PLRun8rBfXkpGLKXoEa7OYHJkNlbALz7/8Z16J6wUm7kCdWLHCY1FU6flj3fWbOkFrSLA/DWb6JZ1PADZl54mB4XxkuvugkTrsiC9lZly0pr3KBmdbFhYGBiJiXBp59KurC3WbVKUtb/vlil7Omn+WFqHurXhyP/mZTotJ7CTXZh+MhqraMrh3kKDcJKuZEryhS6xOrVUK+eVE9IS5OWhHnyuHQIB6ZHU7l3C95MeYfQpN+lJNbIkVlKBPNUFgs0++ABzDwBpGHBhg/MmkXayNHuHpr9TpyQ1ZKVKy/9bkybjZn9onj2Wdk1KPX0CgIrnbjmS68+fjU3No6GI5ZSKXwhDUcs9coXrxqElXK05GTpCrBhg7xvmtCkiVSKeOkl+PhjyTjZutUjqhdly7FjsvfbqJHMggcPBl/X73IdOwbTX4jCnxQs2CQZrFKlHBWALwkJwbIsEuvQYXzSbgXP8D1NJvVk0yZg715pqeiJTBMmT5ZOWd9+K/vZeSQLPMnmz+jYUIYOhT/+gPLlrh+aLj8ullNWkXRPWClHGDcOVqyQjKDdu+UMamgoLFsmgaBqVcn8jIqSjkEAvXpRpsqTHD59nhVfP8eWklVYXbE2qyrWJq2iF5zBWLtWDuBeuCClBd96y+6zrI50/jy0fdikd8I/+Pj7gZUbN1/IKUJCyBMSwuvAokUN6dUL7q1n42DxxyjhdwpjzBgp8bh8uWckbx08KHXCFy2SsTRqhPW2O5j9fCSbx0URWziU92eG0Ly5PNyegiKeUAPbETQIK5VVcXEQfPE/+++/S/CtWVMKzNeocWUbl+++y/h7QgIcOgR58jA4IYD/Tf+T9WXvpMGBrTy0ew0AF0oFQ9XRnlkS6Px5Cba1aklv3TfegDvucMtQ0tKkrHDTDZ/Sy5wAXbrLWDwh8LjIQw/Bli3Qr58PYdM/54fAF7gtLExWA0C2BSIj3fPziI6GMWNg/nx5/7PP4MUXiYm18EIDiIkJoVUr6XxUsmTGl9lTUMTrV5Eu0iCsVGYdOABDh8JPP8Gff0KdOvIk42Pn7k7BgnDnnQCXOvSMDHqLw6cvcK/1FEP846i9KwZKlJBPrlgBL74ILVtCq1ZSODfQdVWcLjl6VALuypWy3J4vH/z4o+vHcZFpQt++UHDhFEYyGDp1gu+/t//3kIMULSpVpKY/2pBGL2xgUvIjtLL9hgGyPRIV5dogbLXKPdu1yyg9OX068Q88xluvwBdfSNCdOlVeRF1v1+BW5T89ocKZI2gQVspeJ0/C8OEwfry8378/lC8vf8/GE79dtYaLF5f7jhoFAQGyx/z991CmTJbva5fjx6VI7+zZsrRpGFL32QOycT/8EPZ9E0mEzzPQJFReEOTCAHy5zp2hSRM/Pun4Nk2io8hDEqZvHixNQzF+/FECYteuznsRd/y47Pd++aWcEU9JAZsN02Lhr5l7aP2iPOTFF6U0dHbaSHtCDWxH0GIdStkjJUVq5R05At27w//+lxGAXeXCBQmEixdLNvLq1XK0Y/hwWQpv0kRm2Hfckb1nt927JdCXLy9nUFu0kLrPYWFSTrDatU9yrj7rPGmS5INtCG5PncL7MVau9M4+dllgz8/aNOGPYdFsHhfFLydCSbknhIV+j1Bq3TwoVkzaD917r6xjZ3Xp/mJVL0JDJRlv3DiYPl3+rzRrJpnyb7+NmZJCsulPM1sktntD+PJLaQvsCJ56xl4rZinlKMePyywUZL2vZk158yTh4dKb99SpjI+FhMAa2V9m/nzIm1eyhZOS4OxZ6WZUo4Z8fsIEmeWfOAERERcLLQ+S4z1paVL/uXr1G94+PUv16hnJdVv4OcCSJbIP2rQpLJqdhP+FeChVyuH38USZ/VmnpUkHwPfeg7//Nnm++gr+FzSK4tHzZKnaMOQFV0SErGnffntGdvvlQfbyIG2a8kIwvR64v79kx0dHywvUvn2hRg327YPf3onmv2lRrM8bStjIEJ57Lkcc2b4lDcJKOUJ0tLya/+gjOXPqyaxWOZ6ycyfs2CFPjAMHyucqV4b9+698fIcOUs0KZO/5+HH5mpAQmfE++iiUK2fXrRuOWHrdvbngoEBWhzfPznd1jU2b4JHGp/jMbxDNYz+lUIXCDr2+p8vqzzo1VWovDxsmicpTg1+lc9xoDEyJiv36wdixEpBr1IDSpSUw22zy7+LTTyWz+d9/5S0hIePiFoscTRsyhAQKMnOm3GvVKonx6b0Y0lMccgMtW6lUdv3+u2Q5lykjCVGezmKRZeJq1WSGcrlVqyQ4//uvzIgLFLgywO7cKUlWWSyu4aos1cOH4YOWUaxKfIaySXEYcb2gQvaaQrhDdpZQs/qz9vOT15Hdu8PEifDz0Mdoz5f4k4LN8Gdd/jbcPa4uBff/BZs3k/bHUnxTUwFIS0rmUMRyKh44IKsp6Ue/vv4arFZMf3/Wl2zP2BcKMmcOJCbKP8MPP5TtZztfy+VaGoSVutrUqfJsVauW9DS9/OyENypT5uYJXEWKZO/yLshSvXAB3m6+mqknHsCCFcPf3ysLcVy9nJxeYALsq4mc3Z91njwy6e3RI4QFb0aSMC+KyYdCifpQlpurV4dyd54nT/1fmLGmN762NFItvryW936av/8Z9UqX4vRp2fnwD3iSPGuj+GpHKPMHhBAUBE8/LXv1997rlb8et9DlaKUut3u3PBM1bizt97KT4JRLOHtP2GaDJ56AfjOb0pQV8kGLRdZWhwzJ9vVdKbtL9874WaekSHG35cvlNNzipWlYk31pQDTNLJEsszZnLddfcbBY5NTcM8/IaaSAgCvH6olJU66gy9FKZdXtt8te6YMPXvn1SOPEAAAgAElEQVSMom7InsIK2fHuuxAx8wzfFdgNib6SGOSlFbGyu3TvjJ91eipASIjk+FV8bTHJRwuy62ARtp3pik+eVAoHbMcnMJVve99NkSJQuLAsoBQpcv3/Jtmd8ecmdgVhwzBKAr+YptnYMAw/YDZQBPjWNM3vbv7VSnk4m02ygdu3lyf29u3dPSKvY9dZ5yyYMkUmvM8+W4h8IzbDrp2yx+2lFbEcsXTvrJ91uuAiAcT5nCFP6TNXfjwo8Jp0gxvJKSUlXeGWJ9sNwygMTALSi8K+BGwwTbMh8JhhGAWcOD6lnCs1VbJHRo+WM7HKY0RHw/BndjG17CC+Gp+GUaK4bBMMGeKVARikwESg35VndDytwIQjxphTSkq6gj3lZaxAZyA9Jz0UmHHx7ysAu9a9lfJIr7wiiVgjRshhSuUR/v0XerQ/wQLzYTol/4T/qf/cPSSH6FAnmOFhtQgOCsRAZpfOOk+dVY4Y441m9h5XUvLcOZg5U16Mu8ktl6NN00wAMDJS3fIB6b2iTgHXpI4ahtEb6A1Q3tVVhZSyV3p5vcGD4fXX3T0addHZsxD2cDLfnn6UcpZD+MyLko72OYSzl5MdIbtj9OiSkgkJklGWL58caH7pJSk+4qajiFkptHoOSH85k/961zBNc4JpmvVM06xXPL3akFKeZu1a+Y83fLi7R6IuslrhqSdN+m9/jvutq/D5cRI0aGDX1+aEBu85hcfN+E+fllqn7dpJBbzp0+XjnTtLWnhzxxaVyYysZEdvABoBvwB3A2sdOiKlXGXCBCnjmBvq6HmJL7pF02zhTJ70nQHvvi9PknbQbFzP4xEz/nPn5PzUr79Kic3y5eWg9L33yueLF88oS+smWQnCk4BFhmE0Bu4E1jl2SEo5UUqKNBd//XUp6+OOloDquqb3j6bn1BbkMVKw+PpJEwA7aTauuoLNJh218uWTzlH9+0ury3r1PK6KiN3L0aZphl7881+gJbAaeMA0TevNvk4pjzJggLQA3LTJ3SNRl5k2DRLGfEcASVhMqyTKLF9u99drNq4C5Az5zz9LN7G4OAm48+ZJ8er69T0uAEMWi3WYpnmYjAxppbzDd99JN/FBg+xe5lTCmdWPIiNhVLeNLPeZgmECPpZMF+NwWYP31FSZYVkssGuXvFDw85Px+vnJW/PmWmnNHTZuhJdflhaf9erBmTMQHOyRgfdyWjFL5Q7r1kkP1Qce0ESsTHLmfmtsLAx6ZA+RZmvylCmGMX4sbN+e6WIcTs3GjY+Xhh7z5kkt8V9/ld7Na9ZAnz7XPn7zZqk7PmWKHH959FFo2zbbNbrVDdhs0j5xwgTZ3/32W9kH9slK3rHraRBWucNHH8mr4mnTMvqlOkBuqI/rrP3Wv/+GZx48zIKkBylUyMTnj8XX7wJlB6eUzoyLk44Ey5dLUk/x4hJQC19sn/jEE1LeNDVV3lJS5M+qVeXz587Bn3/C3Lkycw4Nla/v29fjZ2dexcdHAvGAAfDOO163CqENHFTukJQkvfAqV3bYJV3dzD6rsvtCoVL4Qq73LGEA+0c8nKUxHTsGDRvCg0cmMZaXsEQtlSVEd9u2TXovt20rAbVx44xSpvfdl/lMepsNYmKkHvmcOZIIGBsrn4uNhbvu0uz8rEhKkkLXTz8NderIXrAHvbDRBg5KpZs6Fdq0gaAghwZg8I6MXEcsJTt6v/XsWXjoIZlodlv6NJZKrd3fLjIxUYpUjxwpPXPbtJH93bXZPIHp4yPHYe69V7ZB4uPl40ePQqNGcq8RI+Dhhz0qiNyKW1eA9uyRnI7YWCniUqeOV/3sruYdi+ZKZUVUFDz1FHz8sVMu7w0ZuTd7oWAvR9Y7TkmBxx9No//GbkS+s1zqcLg7AC9eDDVrSpDs2lX2ep01Ow0Kkj9LlJDiESkpUkAiNDT7Ad9F0l/YxcUnYpLxws4lxVGmToV77pG6pvPmSZKll9MgrHKms2fh2Wdlf+7NN51yC2+oj+uIFwqOqn5ks0GPZ006Rz5HV3MyIUE7MvX1TvHXX9IQ19cXli2T42vFijn/voYBjz0my99ffCGZ1vffD//84/x7Z5MjXthlyezZ8qL67rvliGG7ds69n4vocrTKmV59FQ4cgJUr5cC+E3h0fdyLHLWUnN3qRykpsn330LRudONn6NlTiqa4g80mx1nq1ZMn9BkzZM83Tx7Xj8XPT7L2u3WDJUugYkX5+A8/yIuD0qVdP6ZbcPkKUFqavEhq1w4+/xyee05+bjmEzoRVzvPbb/DNN9KY4f77nXYbj6uPex2e0Drv7Flo+5CNdtOepBs/YxqGHN+JjnbZGC45elSWfu+/H/bulY89/rh7AvDl8ueXzGmAI0fk6NMdd8hxGycnz2aWS1eAJk+GGjXg+HEJvH375qgADBqEVU50111SH/Z//3P6rTrUCWZ1eHP2j3iY1eHNPSoAg/tfKBw7JrUrli2DVuV2gGFgmKZMjaOiXDKGS/btk5TsDRvg66+hShXX3t9epUvD1q2ScNSrlxyD8qBlape8sEtOzlghKFVKOnvkUHpESeUsHnZUITfbvx/CHkjgfFw8o34pT9tCK2WJNSVFKkxFRmaqIEe2bNwoGc9WKyxcKMeNPJ3NBl99Ba+9JjPl/fs9pta5U7Oj//1X9stjYuR7/+ADh57tdwU9oqRyp19+kf7AM2ZA0aLuHk2u9tdf0K3lf0w+1Yaq5VPJ23oT+DaWwBsVlemKWNm2YIEEsIgIKQjiDXx8ZPn1oYckESkwUF5kHjkCZcq4dWhO7ZAUHg67d8u56g4dnHMPD6IzYZUzHD0qx0wqVpQjJnbsG+WGalfusHw5DHh4N3OTWlHW/zg+s3+B1q3dM5gzZ6SCkmlKT1lvLx05ZYosUX/4oTSjzymFPmw2+V0VLgwnT8KpU3Dbbe4eVZZlZiase8LK+5mmZNqePStnL+0MwG4765iDzZkDb7Vcxx9JDQkOOo/P8mXuC8CffQa33y71MQ3D+wMwSM3qZs2kRGOTJrBzp7tHlH0nT0qxkrZtJRO6aFGvDsCZpUFYeb/Jk6U+7wcfSAszO7jtrGMOlZwM3z0XzZ9hHzLR73kKlSuIZe0aaR/naqYp/aIHDJCqVG5eunWosmVlaX3SJNixQ45YjR/v7lFl3fr1Unxj6VI5w5ZTZvaZoHvCyrvZbPDpp/Jk27+/3V/mzLOOpik5NOvXS/3+C5HRVPwniu0lQjlcIYTixaUeRLFiXPp7yZKSDFugQLZv73IrVsCEbiuZcOBB8pCKj80PY/wvGY0MXCk1VZZrf/xRsmvHjct5T+yGAd27S5LbK69454uMpCR4910pE1quXEb7wVxIg7Dybj4+EgUSEjL1ZOvIesjx8VITJD3oxsTIChvAg76RzLc+hK+ZijXBjy+TP2DFzsYsSKjK/oSM5LEGRNPciOJ4jVCKtw+heXM5ympPMqy79rZPnJDk1d3fr2KOz2MEkoQBkIq083s4a80dsmXkSAnA770Hb72VszPlS5aUrmDpRo6UpK1hw5xWoMZhbDaYNQt69IBPPvG6zkeOpIlZynvFxsrycxYKLTiiA1JKiqwE/j40mrrnlrGX26lQ3qRF0AYsIfdSrHcYd015Hd9Pr1O7etw4Uvu8SPya7RRu3wjL2XgwIc3w40mmMssMw99fEoibN5dtwPvuk5M9jv4+Mss0Jc69P+AkQ+Jfp4f5LbYSJfGJPy1HgFx9/OhyiYlyBOmxx1x/b3cbMED2wStVknPQLVu6e0RXSkyEUaNknHnzSg6HNy792CEziVkahJV3OnoUqleXIwzffZelS2R1BmmasgU9eDAE74tiKQ/gg5VLcy4/Pyks/+GHkqndvLkknPj5yfJo6dLy4qFSJVm37txZptCXiXn7V6Yntic24hgHtpyhKMdpnWc5KQ1DqfZ0CG3ayFJ2wxFLrzujDw4KZHV48yz9XG5m505Z5U2KimaRb3uCiMcYOFD6uG7enK3jR9ma0U+ZIkd50hsk5FYrVkhZx927ZY/144+lWYS7rVols949e+QI4eOPu3tETqVBWOV8nTpJF5W//nLpuc+NG2HgQFi7PIkqdwYwp/6H3D7pYoMIHx+JUJ9+euXsPDr65sEpOhpatJCpta+vVPsaMkQ2iz/5BAYPvtTPNxV/HmEOEcZD3Hcf7PTbRWDVo/gVP3vFymt2ev1eLS1NnkM3jI/mzJylrM3bnC7vV6f7smcxhr0HtWpl+x7ZmtF//708wb/xhiTn5XZJSRltGdevh9q13VfE5vx5aaAydixUqCBlOJs7/sWhp9EgrHK2uXOlzu4HH8gTrwvExclzyaxJ53gz7yhe8RmH318b8D0alxFAs7MMe6NA/c8/snw3d+6lD9l8LIwYksDcxXk59OdhjlGCBnlX0aLAb6ytUIvN1SpS8bY0ot8OzfL3m5AgdS1+/RUWL0ylc/xXjGIgvqRBQCDG0sx/nzeb6WZ5Rj9/vvxbaNFC/n71en1uduRIRgOInj3lZxMeLsHQVXr3ljruL74orSLz53fdvd1Ig7DKueLjZSm3RAlZwnVyMXebDX58IZp/v4skv+0MvQN/osD5o9Cxo8x4K1S49Uw3u66eKQ8aBO+/D8B/te4j346tBFqTMDBJxZ/W/MYq31Bq1jCoW1eSTuvWlVomFovMbNPfrNaMv1+4ICdF5s2TWs93p6xnsP9YHmIh+VLiMZEZNhaLzLSGDLH7W7jVTLdS+EKu90x00xn9qlWy71mrlgw8lzzBZ5ppSmGPCRPk7927y+/OGdnrmzZJosSgQdKAYuNGOHdOzjTnIhqEVc7199/w5JPSg7VuXafeKjUVPmwXzeCIFgSSKAGodm0pjdmggVPvfY0bBfrZszn91rsE7dhyaU/6aInbGdNzFxs2QO1V49l8oQq+pFKTrUTRjLWkf71JXi6Qh2T8SaExK3iKKSwt1YWAro/zbKHZVP+sN0a7drLk/957WZ7x32qmm+mZsM0mZ7qSkiQYFy9u91hyrYMHZY/4m2/kH/ekSdC1a/avm5oqvX7Hj5ffRWCg5Gk88UT2r+2lNAirnM0F+1tJSbLtXGP+cD4w3sbHtMqe77BhLlsCt9vlM2WLRWboL74oa8qXHf1I/5++qcHzLHv8S4LOHqTHu+WvvZ6fn9SeTC+0kV48Pxsz/lvNdLO0J3zokEzlXbm8mhMcOSL/RgYMgOBg+OknGDNGXliGhMiflStf//+Y1SpJkWfOSGJkWpr8uXevfE2/fvDss1J+MhfTBg4q57lwQQ73h4c7vfzg2bMQ1i6VsOUvc2enWvjM98+YATZr5tR7Z0lIyPUbIxQsKMv3r70G33wjLQQNgzrVk6kzEDhXGPJ+LN/XsmWyDm2aMsuMiro20IaEZHm5/VbnstMD7S2zo48fl5WIN9+U6lEq80qXloS/dAUKyIu1SZPg88/lY8WLS+vHAgUk92LBAnnRc+SIBOLKleXzvr7w8svyfps28kJVZYrOhJV3GDxYnjiWLZNA4yQnT8KjD54nfGMnHmIRjBgh+1nu6PzjKJfPlG+0lGzPY7LBIeeZExMls3bTJqmIUqOGw8ankOC6bRusXStHnNID9dtvy7+PsmXlLThYGqW0aePW4XoyXY5WOUtMjFSq6NlTkkucJC4OOjc/zqg9balvxGB88QX06eO0+7mUPUvJTk4wy9Y5YJsNnnpKzpjOmiUZ0Up5KA3CKudITZX03uPHYft2pxVj2LdPAvC0gw2p5HcQy4xp8MgjTrmXyoJ33pH9+I8+kuV1pTyY7gmrnGPUKKnENHeu0wLw1q1y0iU1uSgFwx7AMrCrFG5WnuHgQVka7dFDtiWUykE0CCvP1rWrJH84aVa6b3I0G3t+y0P5HmLgyjBK1PjCKfdR2VCuHKxbJ8ekcnJDBpUr6XK08kxWqzzhOjHbMiEimjytm+FPMvhYMFat9M7Eq5xq3z5pcde9u7tHolSmZGY5WvPJlWf6+GN44AGpPesENhtEPj8Df5IxuDjBiopyyr1UFsTHQ9u2cpY1vS+kUjmQBmHlebZulTPBxYo5rS/qx+GnqPfPTHnHYpFjOU48+qQyITVVWhHu2yeVmIoWvfXXKOWldE9YeZbUVHjmGSkekF44wMHmzYOgkW9Q2ucYfP4VnD7lvWeAcxrTlGpfkZHSHalpU3ePSCmn0iCsPMtHH8GGDfDLL06pB7xrF3TrBrVrf0yPYWH4tn3Q4fdQ2RAdLWfBhwyRF2NK5XAahJXnSE6W0nlPPCFdihzs7Fn4+ME/KOjXkJ9+LYh/eQ3AHuf++6UqWi7ruqNyLw3CynPkySOzYKv11o/NJNOET9ot56sDbTjy+CuUL//Jrb/IQ2Sr0pS3WLFCMuEbNdK9eZWraBBWniEyUp6ACxZ0yuUnDNnPS8s7klCiKuW/edsp93CGq2sux8UnMmT2FoCcE4i3boX27aUJQEyMNgFQuYr+a1fuFxsLrVvD0KFOufzSX89y/0ftCfC3UWTlvCva+3m6kRG7rmh6AJCYamVkxC43jcjBDh6U332+fFIVTQOwymX0X7xyr5QUePppOY7khJrAR2ZHU7RjU6qzHcusmRi33+bwezjT4eu0/7vZx73K6dPSiefsWfjtNyh/nd7GSuVwuhyt3GvYMNiyRc4NObhPsHVVNIUfb0FxWzIWP198i+Z16PVd4VZ9eL3a559Ly7yICLjrLnePRim30Jmwcp8//4Thw2Um3K6dwy8f8+4CfG0p+GLDsFm9siLW4FbVCPSzXPGxQD8Lg1tVc9OIHOiNN+RIUrNm7h6JUm6jM2HlPgULwoMPwujRDr/03k3nCI6chGEYmD4WDC+tiJWefOXs7GiXZWCbJowcKcfQypeHunUdfw+lvIgGYeV66U1DqlWDRYscfnmrFba2HkR7DpPwwXiCOOPVFbE61Al2aia0SzOwR46E11+HxESnJeIp5U00CCvXGzMG1q+H776DgACHX37+8wvpcPRrtj70GjWH9HX49XOam2VgOzQIjxkjAbhzZ3jbe46JKeVMuiesXGvtWmnMfuGCFOdwsL/XHSdkYk/2F7iLGrPec/j1cyKnZ2CbJrzzDvTvD2FhUhVNjyIpBWgQVq508iR06iRN2r//3uEN2m02+KjP36QYecg3ZzJGgOODfE50o0xrh2VgJyXBwoXQowdMn+6UF19KeStdjlauYbNJc/ajR6VRe+HCDr/FuHEw4a/7aPTtXrq18HPINXNDycjBrapdsScMDsrATk2FtDQIDJR60AUKOPyFl1LeToOwco19++Q4yujRUK+eQy55eYCscTKJu3/YTbs2g+n6rOMCcI4vGYmTMrAvXJBVD9OE+fOdVo5UKW+nQVi5xm23wc6dDmtPeHmANKxWwn/8hLvStrHoqZYYhmOCvMsSljyAQzOw4+Pl3Pfq1fDll7r/q9RN6P8O5VxHj8LYsbIcXaKEw5YjLw+Q3ecup1HKWsLvGsy4QwkOuT7k8JKRznL0qBTfWLcOpk6FPn3cPSKlPJoGYeU8Vit06SLHUv7+26GXTg+E7det4u29n7E84H4Wtqrj0ADp9ISlnMY04bHHpBTl/PlyFEkpdVMahJXzDB0qLQo//xyqVnXopcsEBVLn4A4+i/oICzbuTY2h7pGdDg2QDisZabXCf//BmTOSqJTTpKZKBrRhyKrHkiXQqpW7R6WUV9AgrJxj/Hj44APo2ROefdbhlx/cqhr1Nx7EhoEB+NnSaBS37YoAOTc2joYjllIpfCENRyxlbmxcpu7RoU4ww8NqERwUiAEEBwUyPKyW/XunUVFSlrNwYShdGoKCwM8PduyQz3/3HdxxhySqNW8OAwbAzJmS1OQt/vwT6teHt96S9+vUgfvvd++YlPIimpilHO/gQXj1VWnU/tVXTjmW0qhEMcYd7ER/4xsgmTRfP+p2f5SmFwOkozKb7UpYOnJEOgFFR8OaNZIB/sADMgM+elSW5GvUkLaN589DyZLydcWLS/eg8+fh1Cn5WY0ZIy3+AGbPhkOHJKjdfbcEcE9x9qxUvRo3DkqVgkaN3D0ipbySBmHleOXKwdKlcM894OuEf2JWK/F1mtHyQij//RRJ5QNR+IaG0vSy2tAuyWw+fx7efx8++0yWYwsVkvrU6cGyRQv4668bf327dld2j0pNlQzyQoXk/VmzYMoU+XtgIDRuDI88An3dXIpz1Sp46il5gdC3L3z4oR5BUiqLNAgrx1m5UvY+H38cGjZ02m329h9P1ePRVGrXj8pdQoBrGzO4JLPZ1xdmzJBkpMGDoWbNa47jZKrYh58f1KqV8f7PP8OIETLDXrVK9lpnz84IwqNHw513QtOmTqnBfY3UVBlj8eLyNm2aLj0rlU2Gmd7Rxknq1atnxsTEOPUeygPExkqnorJlYdMmpy2dpuz+h7Q7arAuIJT7ji8gb77rL3U3HLGUuOsE3OCgQFaHN8/azW02OXbz9dey/BwYCOfOQf7813341UviIIldmdpXvtr585Avn+wblywp9w8MlD3lli2hTRu4/fasXftqpgmbN0um8/z5srf9++8Zn9PqV0pdl2EYG0zTtKtggSZmqezbvVuyYQsVkidpZ+1dmiaHHu6N1fTB/PKrGwZgcGBmc7rFi6X3bdeuEviOHJGP3yAAw82XxLMsXz75M29eOHZMWkH27CnJXv37ZwTJuDh48UX48UdZ4rbZMnef0aOhQgWoXTuj41Hzy168aABWyiF0OVplz6FDMgMDWS4tV85ptzqwZBcl965mSp1PeO7pm9/HYaUY4+PhyScluFWqJHu0nTvbVQXK6UvigYEy823TRhKk4uIylqV375ZuRZ9/Lu8XKiRJYF9/DdWry/czenRGMDUMSSSbNUtqPPv4yIuOoUPh4Ycl+Uop5XAahFX2zJwpgWrZMqhm/ywzs40RTBP6jL6D/fl2EflrGbvu4ZBSjPnySdefUaNkLzYTHYDKBAVed0ncacU+gi/7Xps1k9/Lzp3Su3n9eti2LWNGnJoq55ZBfrjpy8v79sns95VX5E0p5VSZ3hM2DKMw8DNQAthgmuZN69LpnnAOZbWC5eJy78GDmZoBZ2WvdNnQKJq/15QxYwxefjlbI7fP+vWyHFuyZJb3P52yJ6yU8njO3hPuBvx88QYFDEdVy1feIzpaMoF37pT3M7kEndm90vM/zabZe80YUmk6/fplacSZM3u2ZBynzwSzuP+Z7WIfSqkcLyvL0SeBmoZhBAHlgIOOHZLyaN98A/36SeDNYgnGTO2Vnj5N2vP92EgdOk7peGny7RSmKWd+X30VGjSQfdZscmh3IqVUjpOVmfAqoALwMrADOOXQESnPlJICzz8PvXvLfuOff8psOAsy0xjhdPvu5L9wlO0PvEzdBk6sGGW1ysx34EDo2FFqXjuo7aJSSt1IVoLwUOB50zTfA3YC1xQGNgyjt2EYMYZhxBw/fjy7Y1SeYPRoyax9/XU5FlOkSJYvZe/xIevoMRRetQAD6LK6ryyDO0tCgmQMDxoE06dL5rFSSjlZVpajCwO1DMNYC9wH/HH1A0zTnABMAEnMytYIlXulpUllqFdekSMubdpk+5L2Hh/aOW0T1QEfTJmJR0VJWUhHSkiQs76FC0NMjJZfVEq5VFaC8HDge2RJOhqY6tARKc9gmvDtt9JQYOVK6QDkgACc7lZ7pQcOQL9NvVhkTMOfFNJ8fFlXugZNHTYCIDlZzsCWLy8lIjUAK6VcLNNB2DTN9UANJ4xFeYq1a2HIkIxWfE4ubXo1c1kUC56NYo3tdZ545GManYphbfla7Njnz/DYOMckOpkm9OolNZmnTcv+9ZRSKgu0WIfKcP68dMeZNw9KlJAG7X374tyU5GvHcP7Jnjx41KBEk7ZsrlaZzVSWzzmyC9KwYTB5snRB6tw5+9dTSqks0CCsZF+0YEGpRwwSmF555aZ1kZ0ledCb5D/6N88XXoDlvmPXfN4hJR+nTJFyjE8/DW+8kf3rKaVUFmkQzs0OH5YZ4ZQp0gCgTBn49Vf3jWf1avy+Gsvn9ON4l+IYPtdm1juk5GOFCnIMacIEbUSglHIr7aKU21y4IEdwuneHKlVg4kTpDOSszkf2stm40PU5/qUCB/uN4H89gh3bBQnkewfpdfzLL+Dvn40BK6VU9ulMOKczTdi6VY741K0LZ8/CE09A0aLQpYssx1au7O5RkpLmQw/jB1JLpDBpRH7yX1wKz3YXpHSnT0sD+h49YPBgB45cKaWyToNwTrRlC2zfLtnNCxdKg4XWreG336QhQWws1Krl2oSrm7lwgZGj8zJ9/70sWJCxFe2wko8pKbL8vG8f3Hdf9q+nlFIOokHYG/37rwSU/fvlz337ZMY7Y4Z8/uWXJQDnzy+9ft9558ozvrVru2XY15WcTHLte+Hvjjz++P94+GEHX9804YUXpNXijz9CkyYOvoFSSmWdBmFPtGIF/PEH/Pdfxlt8POzaJYlEQ4dKw3aQalYVKkCNy45ujxols9xq1TLV/9bloqMxh7xBnj3b2JD3Uz4f44R7jB4N330Hb78N3bo54QZKKZV1GoTd5cwZ2LwZ/vpL3jZtksBbqJD8+cEHcla3ZEkoVQruvFMasfv7Q//+klhVqZJ0M/K96tdYp457vqfMiI6WRhDJyaRhofuLBSld2gn3yZ9f9r7/9z8nXFwppbJHg7A7fP21dCRKV6wY3H03nDolQXjwYJm53Shj2ZOWk7MqIgIzOZn0A0LtC0UBDq4LDdL1qXdvx19XKaUcQIOwK+zeDV99BQ88AA89BG3byvnce+6R4FumzJXnVQsUcN9YXcQsUQKANHzwCfDHp1moY2/w2mvys+3SxbHXVUopB9Jzws6SlgazZ0tiVLVq0iD+r7/kc8HB8NZbEpCDg3NlwYif8vflQSJY1ep9fJZGOrY70rx5MHJkxs9bKaU8lGE6uTh/vXr1zJiYGKfewztg0lAAACAASURBVCOFhsLy5bJn26cP9Owpe7u53a5dnPh1NZWHPUvtOgbLljn4pNSRI9JysWxZaUThyYlpSqkcyTCMDaZp1rPnsboc7UhxcRJoLRapvTxwoMx2r06cyq2SkjA7d8Z32yGKBDzCjz8WdWwAttng2Wfh3DkpxakBWCnl4XQ52lF+/VUKYHz4obz/6KPQvr0G4MsNHozx1190SZvE/8YXpWJFB19/6VKIiJAjWtWrO/jiSinleBqEsys5WY4MdeggR4aefNLdI/JMc+bA+PGM9hlIYMeH6d7dCfd44AE5Y3155rlSSnkwDcLZsXev1CMeM0aqVK1ZA1WruntUnufMGcyePdkaUJfPSgzn668dnIuWlCSlOgEaN86ViW5KKe+kQTg7Tp+WfeA5cyQQ6x7k9RUqxMRGP9AhaRoTfvCnaFEHXz88HOrXhwMHHHxhpZRyLt2wzKzkZJg7Fzp3lif+/fsh0AE9bnOi6GiYN4/1pdrTe357XnwRWrVy8D0iIuQF0EsvQfnyDr64Uko5lx5RygybTfZ8Z8yQLkWa/HNDy3+YR0ivx/CzppKCP4+X+oNp+xqTN68DbxIfL+U8ixSBP//UF0NKKY+QmSNKuhxtL9OEAQMkAI8cqQH4JubGxrH/y+/xs6ZiAD5YaVT2ZxbvinPsjQYNgmPHpJmFBmCllBfS5Wh7ffwxjB0rgfjVV909GrebGxvHyIhdHI5PpExQIINbVbvU+/f7GauYuDUKgDQspPr4sfaOqvwascsx/YFBXhSVKCH7wXXrOuaaSinlYrocbY/t26FmTXjiCZg8GXxy9wLC3Ng4hszeQmKq9dLHAv0sDA+rRYfbg9hZsQbBZ47zknUc5fLvZ2vbQsSWq44B7B/h6IbBSinlWXQ52tHuvBN++w1++CHXB2CAkRG7rgjAAImpVkZG7IK8eYms3oowyy/8XOAJfni6CbHlZOm+TJCDlow/+ggWL3bMtZRSyo00otzM+vUQFSV/b9VKevkqDscnXvMxH5sV499/SDhr8PGR8Sw1H6D4Y39iyZsCyEx5cKtq2b/52rUwZIhkqCullJfTPeEb2b0bHn4YiheHzZu1/ORlygQFEnd5IDZNhi35kra71vD8roc4/E8p3hl3gsXxVg7Hc82ecZYlJ0sjjLJlYcSI7F1LKaU8gEaW6/nvP5n5GobUhNYAfIXBrapdsSf8yuqpdNn0O9PvGMjUZaWYMAGee64YQ2nu2Bt/+KHszy9cCAULOvbaSinlBhpdrpaQAG3awPHjshR9223uHpHHSZ/RjozYRdOo2QxYPYVVNZ7iiW2fMHgwPPecE266e7cE4a5dpTOVUkrlABqErzZuHGzdKrOtenYlt+VKHeoE02H1HFj8BafuCKH5th8ICzOct0pctSp8+aU0ylBKqRxCjyhdzTRh0yaoU8fdI/Fs0dHQogVmUhJJZgD97ohk/IYQx1bESpeSoklxSimvoUeUsiIuDg4elH1gDcA3lpICb78NCxdipqRgmCZ+pPDZo1HOCcB790qLyCVLnHBxpZRyLw3CIDWhn35a2hImJ7t7NJ7r6FHp2fv++5xPsJJk8ycVCz55/CnYLtTx97PZZIP53Dk5q62UUjmM7gkDfPEFREbChAnajvBGYmLg0Ufh5EkOjZxK48+foIJfeyY8FcXtvUMhJMTx95w4UZLjvv4agh1U7lIppTyI7gnv2iXLz82awYIF2hD+en7/XRKiSpVi/RtzefC12gQEyI/LablrcXEy+61bV14g6e9FKeUlMrMnnLtnwmlp0L27dOCZOFGf6C8XHS2z0NBQCYSdOjGl7qc83a841apJ8niFCk68/6xZkJoqqxP6e1FK5VC5OwgnJ0tLwldfhdKl3T0ah7lZhyO7REdD8+by8wkIwPwjkrfL/8gH/aFlS5g5EwoVct74AXj5ZXjkESdHeqWUcq/cHYTz5ZOmDDnI1R2O4uITGTJ7C4B9gXj/fnjpJUhKAsBMTmb6C1F8sDmEXr1k+9zPz2nDhxMn4MgRqFVLA7BSKsfLndnRiYnQuTNs2eLukTjcTTsc3cy5c9C7N9x++6Va2abFQjJ5GLM5lOHDZWXYqQEYpF9zgwZw8qSTb6SUUu6XO2fCb70FM2ZIM4BatS59ONvLuB7geh2Orv743Ng4Fk2YTdXtMfx7+9207NuZDneXhthYeP55CA/nr/kHWPJmFPPPhjJgagidOrlg8IsWSb/md96BokVdcEOllHKv3BeEo6Jg9Gjo2xcefPDSh7O9jOshrulwdNnHQb7P6WOm8f3kIfhbUzFWQK8LiTDoGTqsXcvufRbCX4I5c4IJDg5hxnw5Pu10CQnyAuDOO+GNN1xwQ6WUcr/ctRydlCSz3ypV4OOPr/hUlpdxrzI3No6GI5ZSKXwhDUcsZW5sXLaHnRmDW1Uj0M9yxccu9fLdsoUCT3dl4pQ3CbCm4gOYQO39m/lw1t+81N9CjRpSnOr996VngksCMEB4OBw6BN9+q2e1lVK5Ru6aCU+YAH//DX/8IUlZl7FnGfdWPGE23aFOMD7JSawfO4kGscsobKbi07cvIXWCYfNJ7ty/hfXlatLwn7/wMW2kWnyZn9iBmJEN2GiVbeGhQ6FkSZcMV5imLD+/+qrsByulVC6Ru4Jwnz5SealFi2s+datlXHvcbDadHoTt2ne+/Izu1ZWokpOl3/HixbBsGRQrJqUk27eXmX6FCrQ/doz2l3/NwcbyZ61aPPbGdOLOJFHn0A7uifmXiIMdWb0plMLVj7N6VnGqV7f723Ucw4Bhw9xwY6WUcq/cE4RTUmSZs2PH63766kb1cNkyrp1uNZu2a6a8dKn0y01JAYtFDubef78kk4G8iLg6c3jPHgnCAQGS9b1jh1SZMk25RpEiAFxINGgaWJsvZlxg/p4WzD0fgH+peMqFrWfswOBLAdilCWqjR0ONGlfszyulVG6RO/aEV6yQfrQ3OZLUoU4ww8NqERwUiAEEBwUyPKxWpoLPjWbN6R9PnynfE7eDvtEzuOfQdoKP7Gf5uJ8zHvzUUzLbNU2p6LVmDRw4kPH5YcOkhrPPxV+dxQKNG2d8fuxYeO89CcgWCzY/f349E0q7drLiO3xAUVL2BBNU+QzFH9lA3Zc2MnZg8BUz9SGztxAXn4hJxgsFp+xtb9wIgwfD9OmOv7ZSSnmBnF87Oi1Nyi7Gx8sM0Sn99sTVM12Q2XR6MK8UvpCQ/ZuY9MtQfG3yGAM4FViQIufjZVl2xAjZlLVapYduZOS1S9IXe/le6rN78TEJCVIKe+dOOLckGtvSKCbHhbKWECpUkMly+/bQpMmN2/M2HLH0usvywUGBrA5v7qgfFVy4IL+XhATYuhUKF3bctZVSyo20dvTlvvlGik/MnOnUAAwZS8o3WsotExRI/9VT8LsYgG3AotsbMvnhXkxNv0h4ODRtesM9YasVDpQK4diHkVgjo1jjH8pvb4WwcyccPpzxOIslhLp1Q2j/AkxoDzVr2leC2REJanZ59VV5xbBkiQZgpVSulbOD8MmTspfarNkN94IdrUOd4CuXsLdsgWfehB49GNyqGpM2dOTuOXuw2KykWnyZfH9HOndreUWENBuE8HeJEPbsgb3jpa/93r2y9bt/v/Q1gBAghEKF4I47ZOv4jjsy3ipXvvFs92YckaB2S6tXw1dfyVL0dZLklFIqt8jZQXjSJDhzBsaMcW0nnjVr5Lzrtm2wbp3MwBs2pMNzTeCNPrxUrDhVt8ew9856dO4dRoc6wRw4IKvKf/whb8eOZVwuXz647Ta46y4IC5Pt7Sr/b+/O46Mqrz+Ofx5iWGQVwbBKwSoVq4gVa0QlVYpWrGhrxIVa5CcidUEroIiKoiAVteKC1qVUrFXr8sNaf5W6NKKoEJCtILIjq1JtRBbZ8vz+OBMIMQmTyZ25cyff9+uVV5Lhzs25c8mcPPc+zzmHWbLNyQn20IKYoLZfJ51klbHy84Pbp4hIBGX2PWHvYc4c6xecKtOm2U3X4mL7fuBAGDNmzwzlEv/9r11xLkm6ixfb4zk5tuKoe3dr8PT97wefaPcnabOji4tteVWrVtXfl4hImtI9Ye9tKJmTk9oEDDYTu0RWlnUCKpWAZ860eVdvvGE5qX59u/U7aJAl36OO2n/CTfYSou9cUg/KAw/AHXfYrOjDDgt+/yIiEZOZSbikOcP778Oxxyb/5+3aZU0HfvADy6h16uyduZyXB9jcsJEjYfJky8k33QRnngk//nHV7t2mQ1WuhMyZA8OHQ69edsNaREQyMAlv2QJDhlhLvlIdkpLmiy/goousyMbgwXDppXZzNza7edFBudx+oS2FbdTIlvAOHmxfJyKeqlxpZ+tWe42aNbPZ6qm8ti4iksYyLwmPHWuNAJ5/3i4HJ9MHH8AFF9gs7IkToV8/ezw3l2WH5DJqlM0/qlcPRoywVTnVXY2TsiVEQbrhBlu8/OabalEoIlJKZiXhlSth3Di45BLo1i05P6OkrvPhh1t1q7Zt7bHYZe9vv7WB+B/+AAccAL/9LQwbBs2bB/PjU7KEKEi7d9ul+SFD7Ka3iIjskVlJ+I03bPQ7dmxy9l+2UtXQoZZcYsPbL76Ac8+1zX7zG1ui3LJlsCGkZAlRkLKybLlWyWxxERHZI7NqR195pbUqbNMmOft/5x0b6paM7ho02JOAFy60SVazZ1txrkceCT4BQzA1rlNi+3YYMMBmpMHeWtciIrJHZoyEvbeSUocfnrxGuN7D3Ln2uVatfWY+v/kmnH++1eSYOhW6dk1OCCWStoQoKMXF0L8//OUvVq3smGPCjkhEJC1lxvBkyhTo2NF67CbL/ffbELdPH7jrrj1NEx57DH72M1sOPH168hNwJIwYYQl4zBi7by4iIuWK/kh4926b+dS+/Z6RaeCef97u/ebnW3KpVYvdu2Hob60d7lln2SYNG1b/R6W0l28yPPaY3ZMfONAWQ4uISIWin4SfecaaJLzwQmIdC+Lx+utw8skwaRLUqsXmzTbAe+01uPZauO8+mwldXZEtxFGiuNjOQ69e8PDDWg8sIrIf0a4dvW2b3Qdu1cquBSfrTb+42IqANGzIxo3Qs6fNN3rwQbjqquB+TMp6+SZTycS1+vXDjkREJBRVqR2d8D1h59wE59zPE31+IAoL4auvbG1w0Al43Tq72fvZZzYRq2FDduywCViLFtngOMgEDBEtxAE2Iz0/37pS1K2rBCwiEqeELqI6504BWnjvXws4nqo59VRYvTr4KkybNlkCXr7ckvyhhwJ26XnqVLstfOaZwf5IiGAhDrBqYWedZYukP/+8+iXBRERqkCqPhJ1z2cATwErnXO/gQ4rTwoW2XCjoBDx1Khx3nPUCfvnlPZWwHn3UqmC17r6Sm+e+Trex7zB59tpAf/TQMzpSL3vfUptpXYjj22+hd29YsQJefdUaWIiISNwSuRx9KbAQuAc4wTl3TdkNnHNXOOdmOudmbty4sboxftfy5ZYcx40Ldr8lFbGWLbNKT7Hpzv/6F1xzraf+4V+QdcICPHsnTQWZiCNTiAOsWEnfvtY/edIkOOWUsCMSEYmcRC5HdwEe995vcM79GRgNPFR6A+/948DjYBOzqh1lWSNG2HTkvn2D3e9bb1lbQrDJRQUFrGiRS34+1G66laa9ZuNK/dmSjO5FaV+Io0RREcyYYeun+/QJOxoRkUhKZCS8FChpCHs8sCq4cOJQWGiLcm+4wWZFB6lHD2t5lJUFtWuzpWse55xjk6ObnjuDWnV2fecpaT9pKmgLFtgfKoccYkvDrr8+7IhERCIrkST8FPAT59xU4DfAvcGGVAnvrWlC8+b2OUgTJ8KRR1olrDvvpPjNt+n7SC6ffAJ//Su0a1/+gD6tJ00F7c9/huOPt0pYAI0bhxuPiEjEVflytPf+GyA/CbHs37p1sHgxjBwJjRpV+ekVVqP6+9+t1vGGDTB8OOTmcvttMHkyPPCADZA3Hxyx7kVB2rXLql/dd59VJRs0KOyIREQyQvSKdWzZYpWxsrOr9LSy1ajAkui9p7Wh18U/tcurhYVQuzYvvggXXGB5+ckn9y5BjnxJyUR89RVceKF1qbj6arsHXMXXXkSkJqlKsY7ola1MsBDEuCmf7pOAwSZW1bnuGlvr+sYbULs2c+dCv35w0kkwYcK+NUAiM2kqSKtW2QSsJ5+E//mfsKMREcko0UvCCSpvAlWvT96jx7wCGD0aOndm50741a+gSRN45RWoUyf1caaFoiJ46SVLul26wMqV9qKIiEigMqOVYRzKm0D1cesf8Hy3X1oXJmzZ8fz5VpgjWW2J09rGjbb8q107GDDACmSDErCISJLUmCS8TzUq78F7ig5uQd2HxsMBB7B4MYwaBb/8JZxzTrixptzXX9uSr+99D+6+22pyzpkDnTuHHZmISEarMZejS+7ljpvyKSdN/Rv5S97ni6ee4ewurfHe2t/WrQsPPbSfHWWSHTv2TnJ7/nnrTjF8uMpPioikSI1JwhCbWLViOox8FDp1gpNtedHEiVBQYLWhW7YMN8ak2rTJWj5OmwYffGDNL+bPhwMPhCVL7LOIiKRMjbkcDVjiyc+3EeCiRTB9Op9/DkOGWOnjyy8PO8AAeW+NFXbutO9/9zu7t9uzJ9x5p93/Pfts68kMSsAiIiGoUSNhHn3UalCCJaeCAgaPz2XLFnj8cWsbnPa8tyVV69bB+vXWyCInB2bNspllGzbYx/r1NvKdMQO6doXcXCty0q0bnHBCQsVOREQkWDUrCc+bZwt/a9WC2rWZlp3HCy/AHXek0W1Q72HpUkuqa9ZYua5jj4WPP7ZZY+vW2Ui+xEsv2eObN9s2LVrAMcfYiPfII/f0QubUU+1DRETSRs1Kwu+/byUqV65k6wl5XNw/l06drCJj6P7zH7j2Wnj3XUu0JcaPtyTcrJmNYlu3tsYVJR+dOtl23btbSU8REYmMmpGEN22y7kgNG8JFFwFw83Xw2WeWl2vXTnE8a9bAa6/ZbLCjj4ZbbrFmCDNn2mg1L89KdrVrt6enMYceag0UREQkY9SMJHz99XZ5t7AQsrOZMQMefND6EHTrlsI4Vq60dbgTJ9o96TZt9l4Hz87WSFZEpIbJ/CQ8a5YlvRtugOxsdu60YlAtW1o+TKnhw60e5oABMHgwHH74vsWpRUSkRonCfODEeW/Jrlkzu+SLdeObNw8efjgF7XCXLrVWTAsW2Pdjx8KyZfDII3DEEUrAIiI1XEaNhMu2GnyAT+k6bZqtP2rcmA0bbIls795w3nlJDGTJErjrLnj2WbvMnJcHRx1l93hFRERiMiYJl+0XvLZoG9+88jRFHY+iSf/+gNWG3rHDltMmzaRJcMUVtgzq2mth6NAML8MlIiKJypgkXF6/4MvPvZmjs7bxalYWixfbgHjgQLsVmzQLF9rM5mefVfIVEZFKZUwSLt0v+OAtRezMOoBNdRswzzcA7JZw3bpw661J+OFffmnLjjp3tt7ExcV2GVpERKQSGTMxq3S/4FvfeYI3/ng1tXftpFWTesyYAS++aBOkW7QI+AfPn29lIXv3tmvdWVlKwCIiEpeMScIl/YIvnP0Pzl34Lu+360xWvboM6dmRYcOgeXNr1BCol1+2mszffmutAFNe9UNERKIsY5LwuV1a81iHbxnz5gQ8cM6i93jssB3U/bw1774Lt922t/hUtRUX2w7PP98qXs2cCSeeGNDORUSkpsiYJAzQ/Z8vUMt7HFDX7+aUtQu48Ubo0MEmLAfGe2uW0L+/lZ5s1SrAnYuISE2RMROzACsFWapL0pTtecyfD889F9CV4t27YetWG1K//LLtVAU3REQkQZmVhF96Cd56CwoL2Z6bx5W/zuVHP4ILLgho/zfcAG+/DR9+CA0aBLRTERGpqTIjCW/bZk3sO3Sw/rs9evDI/dYl6Y9/tIFxtY0fbx/XX68ELCIigciMe8IPP2zdiJYvB6CoyJbr9uwJp58ewP5ffdWS73nnJbncloiI1CTRT8Jff22NEU4/3UbCwD33wFdf2cPVVlhoPYi7drV+vllZAexUREQkE5Lwvfdaxh0zBoC1a+GBB+CSS6BLlwD237KlDan/9jc48MAAdigiImKifU/488/h97+3mVexjHv77bBrl3VLqpZvvrGk26YNTJ5c7VBFRETKivZIeOpUWzYUy7iLF9tErEGDoH37aux3+3Y45xwbTnsfTKwiIiJlRDsJ5+fb9ecjjgCshW+dOnDzzdXYp/dw+eVWhOPnP9c6YBERSZroJuGlS+1z06YALFli3QMHDYKcnGrs9557bALWnXfaSFhERCRJopmEFyyAjh3hiSf2PDR6tBWwGjq0GvudN896Hebnw4gR1Y9TRESkEtFMwrfcAvXr27pdYNkyG7xeeWU1WxVu2mQ9gSdM0GVoERFJuujNjp4+3WYrjxoFzZoBtjopOxuGDavmvk8+GWbMUAIWEZGUiNZI+MMP4eKLoUkTuO46AFasgEmTrEtSy5YJ7vfjj+0y9I4dSsAiIpIy0UnCH34Ip51mpSm3bIF//xuwUXBWFtx4Y4L73b4d+vWDp56CzZsDC1dERGR/opOECwqsVSFAcTEUFLByJfzpTzBgQDVa+o4aBfPn2ySv2ExrERGRVIhOEs7Ls+nPWVn2OS+Pu++2DkkJj4ILC63A9GWXQa9eQUYrIiKyX9GZmJWba718CwogL4/PWucycaLV1WjTJoH9eb93CH3//UFHKyIisl/RScJgiTg3F4C7B9lDN92U4L6cg6eftmVJTZoEE5+IiEgVRCsJx6xebfOo+veHQw9NYAebNkGjRrYmWEREJCTRuSdcSkmf4OHDE3jy1q3WG7haBaZFRESqL3JJeM0aePJJW1XUrl0CO7j1Vmu31KNH0KGJiIhUSWSS8OTZa+k29h069V7Jzl3FHN97Q9V3snAhjB8PAwfammMREZEQRSIJT569luGvzGfV6mK+mduW+j9cw+8/msPk2WurtqMhQ6BBA+t5KCIiErJIJOFxUz5l287dbFvSAoodjXKXsm3nbsZN+TT+naxeDdOm2eXoWM1pERGRMEVidvS6om0ANDxuFXXbbyS7ybZ9Ho9L27bWg7hRo2SEKCIiUmWRGAm3alJvz9fZB20t9/FKLV1qpS6bN4c6dYIOT0REJCGRSMJDz+hIveysfR6rl53F0DM67v/JRUVW4GPw4CRFJyIikphIXI4+t0trwO4NryvaRqsm9Rh6Rsc9j1dq9Gj48kur7CEiIpJGIpGEwRJxXEm3tGXLbElSv37QpUtS4hIREUlUJC5HJ2zYMOu4pCVJIiKShjI3CRcVwZw51ucw4WbDIiIiyROZy9FV1qSJVcjyPuxIREREypWZI+G5c61RQ506ULdu2NGIiIiUK/OS8JYtcNZZ0Ldv2JGIiIhUKvOS8D33wLp1VidaREQkjWVWEl6/HsaNgz594KSTwo5GRESkUpmVhMeMgZ077bOIiEiaSzgJO+dynHOzgwymWoqLbTb0ZZdBhw5hRyMiIrJf1VmidC8QZweFFKhVC956C7ZvDzsSERGRuCQ0EnbOnQZsATYEG06C1q+HDRvAOS1JEhGRyKhyEnbO1QZuBW4KPpwEjRgBRx1la4NFREQiIpGR8E3ABO99UUUbOOeucM7NdM7N3LhxY+LRxWPJEpg0CS69FA48MLk/S0REJECJJOEewFXOuQLgWOfck2U38N4/7r0/3nt/fPPmzasbY+VGjbImDTelz8BcREQkHlWemOW9P7Xka+dcgff+8mBDqoJPPoFnn7XCHDk5oYUhIiKSiGqtE/be5wUUR2LefhsaNrSWhSIiIhET7WIdV18NK1ZAs2ZhRyIiIlJl0U3Cq1bZ56ZNw41DREQkQdFMwrNmQfv28OKLYUciIiKSsGgm4ZEjoUkT6Nkz7EhEREQSFr0k/NFH8PrrMHQoNG4cdjQiIiIJi14Svu02m4h1zTVhRyIiIlIt0UrCq1fDtGlw443QoEHY0YiIiFRLdboopV7btrB8ua0NFhERibhoJWFQZSwREckY0bocLSIikkGUhEVEREKiJCwiIhISJWEREZGQKAmLiIiERElYREQkJErCIiIiIVESFhERCYmSsIiISEiUhEVEREKiJCwiIhISJWEREZGQKAmLiIiExHnvk/sDnNsIrApwl82A/wS4vzDpWNJTphxLphwH6FjSUaYcBwR/LO28983j2TDpSThozrmZ3vvjw44jCDqW9JQpx5IpxwE6lnSUKccB4R6LLkeLiIiERElYREQkJFFMwo+HHUCAdCzpKVOOJVOOA3Qs6ShTjgNCPJbI3RMWERHJFFEcCYuIiGSEtE3CzrmnnHMfOuduqc42YXPONXbO/cM590/n3P8652qXs80BzrnPnHMFsY+jw4i1MvHG6Jy7wzlX6Jx7JNUxxss5N6jUccxxzv2hnG2icE5ynHPvxb7Ods695pyb5pzrX8lz4tou1cocy6Gx1/wd59zjzjlXwXNaO+fWlDpHcS0JSbYyxxJ3jOn2flbmOO4odQyLnHPDK3hO2p2T8t6D432tU3FO0jIJO+d+AWR573OBDs65wxPZJk1cAtzvve8JbADOLGebY4DnvPd5sY/5KY0wPvuN0Tn3I+Bk4ATgC+dcj1QHGQ/v/aMlxwG8BzxRzmZpfU6ccwcBTwP1Yw9dA8zy3ncDznfONazgqfFulzLlHMtAYJD3/jSgLVDRH0A/BkaXOkcbkx9t5co5lrhiTLf3s7LH4b0fWep35t/ApAqemnbnhO++B19IHK91qs5JWiZhIA/4a+zrf2Jv7IlsEzrv/QTv/Zuxb5sDX5Sz2YnA2c65GbG/vA5IXYRxiyfG7sDL3iYaTAFOSWmEVeScaw3keO9nlvPP6X5OdgN9gE2x7/PY+/swFahozWO826XSPsfivR/hvf8k9m8HU3ERhROBy51zHzvnn0IM6wAAAp9JREFUxiQ/zLiUPS/xxphHer2flT0OAJxzXYE13vu1FTwv7c5JOe/BfYnvtc6Lc7tqSdckXB8oOclfATkJbpM2nHO5wEHe+4/K+edCoIf3/gQgGzgrpcHFJ54YI3VOgKuARyv4t7Q+J977Td77r0s9FO9rn3bnqJxjAcA51wdY4L1fV8FT/4G9UXYFcp1zxyQvyviUcyzxxphW56WicwIMBh6q5Klpd05KlLwHA6tJo9+VdE3Cm4F6sa8bUH6c8WyTFpxzTbH/uBXdg5vnvV8f+3omkI6X1uOJMUrnpBbwE6Cggk2icE5Ki/e1j8Q5cs51AIYA11Wy2Qfe+2+897uB2aTnOYo3xrQ/L865JsAh3vtllWyWluekzHtwWv2upN2JjpnF3qF/Z2BlgtuELjYR60VguPe+ohrazzjnOjvnsoBzgbkpCzB+8cQYiXMScwow3Ve8Ri8K56S0eF/7tD9HsfuRzwH9KxiNlZjinGvpnDsQ6Indq0w38caY9ucF6A383362SbtzUs57cHr9rnjv0+4DaIS96d0PfBJ7Ae7azzaNw467gmMZBPwXG3EVACPLOZYfAvOA+dikhtDjLuc49okRaAo8WWabWsA0YDzwKdA+7LgrOZ4xwC9iX3eK4jmJxVkQ+9wOWBB77QuBLOA04Ooy239nu7CPoZxj+R2wvtTvTPcKjuUnwKLYebo61fHGeSzfibGC/29p+X5Wchyxr/8CHFfq+0ick3Leg39d9rUO85ykbbGO2F/DPwWmeu83JLqNpJZzrh7QC/jYe7887HhqEudcK+wv9ym+khFkvNtJaun9LHXifa1TcU7SNgmLiIhkunS9JywiIpLxlIRFRERCoiQsIiISEiVhERGRkCgJi4iIhERJWEREJCT/D0B0cJDjGaPuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(x, y, 'o', label=\"data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "ax.plot(x, res.fittedvalues, 'r--.', label=\"OLS\")\n",
    "ax.plot(x, iv_u, 'r--')\n",
    "ax.plot(x, iv_l, 'r--')\n",
    "ax.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS with dummy variables\n",
    "\n",
    "We generate some artificial data. There are 3 groups which will be modelled using dummy variables. Group 0 is the omitted/benchmark category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "groups = np.zeros(nsample, int)\n",
    "groups[20:40] = 1\n",
    "groups[40:] = 2\n",
    "#dummy = (groups[:,None] == np.unique(groups)).astype(float)\n",
    "\n",
    "dummy = sm.categorical(groups, drop=True)\n",
    "x = np.linspace(0, 20, nsample)\n",
    "# drop reference category\n",
    "X = np.column_stack((x, dummy[:,1:]))\n",
    "X = sm.add_constant(X, prepend=False)\n",
    "\n",
    "beta = [1., 3, -3, 10]\n",
    "y_true = np.dot(X, beta)\n",
    "e = np.random.normal(size=nsample)\n",
    "y = y_true + e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         0.         0.         1.        ]\n",
      " [0.40816327 0.         0.         1.        ]\n",
      " [0.81632653 0.         0.         1.        ]\n",
      " [1.2244898  0.         0.         1.        ]\n",
      " [1.63265306 0.         0.         1.        ]]\n",
      "[ 9.28223335 10.50481865 11.84389206 10.38508408 12.37941998]\n",
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 2 2 2 2 2 2 2 2 2 2]\n",
      "[[1. 0. 0.]\n",
      " [1. 0. 0.]\n",
      " [1. 0. 0.]\n",
      " [1. 0. 0.]\n",
      " [1. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "print(X[:5,:])\n",
    "print(y[:5])\n",
    "print(groups)\n",
    "print(dummy[:5,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.978\n",
      "Model:                            OLS   Adj. R-squared:                  0.976\n",
      "Method:                 Least Squares   F-statistic:                     671.7\n",
      "Date:                Tue, 24 Dec 2019   Prob (F-statistic):           5.69e-38\n",
      "Time:                        14:59:00   Log-Likelihood:                -64.643\n",
      "No. Observations:                  50   AIC:                             137.3\n",
      "Df Residuals:                      46   BIC:                             144.9\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.9999      0.060     16.689      0.000       0.879       1.121\n",
      "x2             2.8909      0.569      5.081      0.000       1.746       4.036\n",
      "x3            -3.2232      0.927     -3.477      0.001      -5.089      -1.357\n",
      "const         10.1031      0.310     32.573      0.000       9.479      10.727\n",
      "==============================================================================\n",
      "Omnibus:                        2.831   Durbin-Watson:                   1.998\n",
      "Prob(Omnibus):                  0.243   Jarque-Bera (JB):                1.927\n",
      "Skew:                          -0.279   Prob(JB):                        0.382\n",
      "Kurtosis:                       2.217   Cond. No.                         96.3\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "res2 = sm.OLS(y, X).fit()\n",
    "print(res2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare the true relationship to OLS predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlcVFUbwPHfHXZXEjUVcMsN963SLCMttdQy23wzzdLMNy1tMTUtLfcss01Nsyxtea2Uyg0VRU0xFzBxQy03cEdxZZ057x9HRIxlgBnuIM/38+EDc2fm3iOD88w55znPMZRSCCGEEKLwWcxugBBCCFFcSRAWQgghTCJBWAghhDCJBGEhhBDCJBKEhRBCCJNIEBZCCCFMIkFYCCGEMIkEYSGEEMIkEoSFEEIIk0gQFkIIIUzi7uwLlC9fXlWvXt3ZlxFCCCFcxrZt284opSrk9jinB+Hq1auzdetWZ19GCCGEcBmGYRy253EyHC2EEEKYRIKwEEIIYRIJwkIIIYRJnD4nnJXU1FRiY2NJSkoy4/JFnre3NwEBAXh4eJjdFCGEEAVgShCOjY2ldOnSVK9eHcMwzGhCkaWUIj4+ntjYWGrUqGF2c4QQQhSAKcPRSUlJ+Pn5SQDOB8Mw8PPzk1EEIYS4CZg2JywBOP/kdyeEEDeHIpGYFRIVR5tJq6kxfAltJq0mJCquwOfs06cPTZs2pWXLlsyePTvbx4WHh3Po0KECX08IIYS4kSlzwnkREhXHiIXRJKZaAYhLSGTEwmgAujXzL9C5P/vsM4KCgmjSpAl33nknjRs3/tdjwsPDCQ4ORqp+CSGEcDSX7wlPCY25FoDTJaZamRIa45Dz+/n50blzZ0JCQrj77ru55557GDlyJADPPfccc+fOZciQIfTs2ROAY8eO/etxQgghRH64fBA+lpCYp+P5kZ4kNmnSJJYtW8bvv/8OwNdff02fPn2YNm0a3333HQBxcXH/epwQQgiRHy4/HF3F14e4LAJuFV8fh13j7NmzBAYGMmnSJEqVKsXFixezfay7uzvvvvturo8TQginUQp27YL69cHi8n0pkQOXf/WGdqyLj4dbpmM+Hm4M7VjXIedPSEhg2bJlbN++nREjRvDll19myj728fHhypUrgF6jO3Xq1CwfJ4QQhSIuDtq1g0aNYOFCs1sjCsjle8LpyVdTQmM4lpBIFV8fhnasW+CkLICXX34ZLy8vJk+ejGEYDBgwgAoVKlCiRAni4uLw9/fnscce4/nnn+e9995j/vz5dOnSJcvHCSGEUykFc+fCq6/CpUtw551w++1mt0oUkKGUcuoFWrZsqW7cynDPnj0EBQU59bo3O/kdClHMTJ4Mw4dD27YwZw7UqmV2i0QODMPYppRqmdvjXL4nLIQQxZZScO4clCsHzz8PZctC//56HjgqCi5cgHvvNbuVogAkCAshhCs6dAheeAEuXoQNG6BCBRgwIOP+0aPh8GH46y/TmigKzuUTs4QQolix2WD6dGjYEDZtgueeyzoDOjAQYmMLv33CoaQnLIQQruLECejRA9auhQcegNmzoVq1rB8bEABnz8KVK1CiROG2UziM9ISFEMJVlC2rg+qcORAamn0ABt0TBukNF3F2BWHDMMoZhvGAYRjlnd2gwvLZZ58RHByMj48PwcHBLFq0yOwmCSGKo5gYePppuHwZfHzgzz91ElZudQgCAvR3CcJFWq5B2DCMW4DFwB3AGsMwKhiGMccwjAjDMEY5vYVOMmjQIMLDw/H39yc8PJxHH33U7CYJIYqTtDR4/31o0gSWL9cVsCD34JuueXNYtw5a5roKRrgwe+aEGwOvKaU2XQ3I7QA3pVRrwzC+MgyjtlJqf34bMGQIbN+e32dnrWlTmDYt788LDg7m9ttvZ8eOHYSGhjJmzBiCg4MJDg5m7ty5ADz55JP07t2bU6dO0ahRIz7//HPHNl4IcfPbtUsnXG3ZAo8+qhOxKlXK2znKlIF77nFO+0ShybUnrJRaezUAt0X3hjsCC67evQK424ntK1SbNm2idevWhIaGZvuYWbNm0bBhQ9atW8fx48fZsWNHIbZQCHFTGDwYDh6E//0Pfvkl7wE4XUgIrFjh2LaJQmVXdrShiyQ/BZwDFBB39a6zQPMsHt8f6A9QtWrVHM+dnx6rszRs2JDu3btneV9iYiI+Pj7ExMSwceNGwsPDSUhIIC4uLst9iIUQIpOoKKhcWQfcOXN0RnOFCgU753vvQZUq0KGDY9ooCp1diVlKGwjsAO4C0rcwKpXVOZRSs5RSLZVSLSsU9I+sEJUqVSrTbU9PT06fPg3A8uXLAahbty5DhgwhPDyccePG5fohQwhRzCUnw6hRus7zqKtpNNWqZRuAQ6LiaDNpNTWGL6HNpNWERMVl+ThAZ0gfPeqERovCYk9i1jDDMHpfvekLTCJjCLoJcMg5TTPfww8/zKeffsqAAQPw8/MD4IUXXmDZsmW0bduWmTNnEpi+TEAIIW60ebNOoBo/Hp55Ridi5SAkKo4RC6OJS0jUQ44JiYxYGJ19IA4IkOzoIs6e4ehZwALDMPoBO4EQYJ1hGFWAB4FWTmyf0x04cODaz+Hh4ZnuS5/7vdGCBQv+dUwIITL58Ufo2VMPFy9ZAg89lOtTpoTGkJhqzXQsMdXKlNCYrHeOCwyUgh1FnD2JWeeUUg8opdoqpV5SSp0HgoFNwH1XbwshhABISdHf27eHV16BnTvtCsAAxxIS83Rc1goXffmqmHU1MC9QSp1wdIOEEKJIunQJXn5Z72pkteo5348+0lWw7FTF1ydPx+naVWdZ16yZnxYLFyBlK4UQoqBWrYJGjeDzz+GOOyA1NV+nGdqxLj4ebpmO+Xi4MbRj3ayfULYsVK8O7rINQFElQVgIIfLr4kXo109vtuDpqStYffwxeHvn63TdmvkzsXsj/H19MAB/Xx8mdm+U9Xww6P2Gp0zRdaZFkSQfn4QQIr/c3SEiAt58E8aM0bWfC6hbM//sg+6NDEMH4e7doWPHAl9bFL5i2xMeM2YMQUFBtG3blvbt23Ps2LF8neP6jOohQ4bk+pzt27ezPQ91Om+8hhDCZGfOwKuv6jlgHx+IjITJkx0SgPMlIEDWChdhRScIR0TAxIn6u4OMHDmSdevW8dxzz/Hpp58W+HzT7Cj/ldcgLIRwEUrBggVQvz589hn88Yc+7uVlbrsCAyU7ughzjeHo4OB/H3vySXjpJb3+rU0b2LEDbDawWKBxY117tU8f/an08cczPzePPcdz585d29Lw+g0crly58q/NGs6dO8cTTzyB1WpFKUXwdW0PDg6+1mtNSkqiT58+xMbG4uvry4IFCxg7duy1LRPnzZtHWFhYnq8hhDDB8eP6/SgkRO9aFBamE7FcQUBAxgcCUeQUjZ7w+fM6AIP+ft4xS5PHjx9P27Zt2bRpE4MHD/7XBg5ZbdYwa9YsunTpwpo1a/Dw8Mj23LNmzaJJkyb88ccfPPbYY+zcuZOJEycyfPhwhg8fTlhYWIGvIYQoJAMG6O0G339fj8a5QgBOHx202eDcOd1hEUWOa/SEc+q5ligB332nF76npOgMxO++g9at9f3ly+e555tu5MiRPPPMM9du37iBQ1abNRw8eJCnnnoKgJY57OO5d+9eHnvsMQD69OmT7eMKcg0hhBMdPqyHmitV0ut909KgTh2zWwVA5IAvaPLFQCwWheHpqd8DpWJWkVQ0esKtW+vhn7Fj9ff0AOxgN27gkNVmDVWrVmXX1c23c5rbrVevHlu2bAFgwoQJfPnllwD4+Phw5eonVqVUga4hhHACm02v923YEF57TR+rWdMlAvDxmAssqzmQ5l8MwA0rhs2m1yRv2GB204quxYth/37TLl80gjDowDtihNMCcFay2qyhf//+/PLLLwQHB3PhwoUcnxsZGUlwcDCRkZH06tULgAceeICFCxfSpk0b1q9fX6BrCCEcbN8+naMyaJB+r5kwwewWATonbOmQFViDGtDx4Az2NHqSRHywGm7g4aEztFevNruZRY/Npjt3kyeb1gRDKeXUC7Rs2VJt3bo107E9e/YQFBTk1Ove7OR3KISDLVmikzy9vGDqVHjuOb0O12QHDsALL0CZ8F+ZVmIk7nO/JPCJVvSoFkHXMuH0/LQ13HefDibpWyWK7CkFP/+sf2fly0NcnP7u4Cx3wzC2KaVynU8sOj1hIYRwBuvVXYvuvBN69IDdu+H5580NwBERWMdN4I/7RzMz6GOioqDLrEeodm47gU/ojese52fKH9+pe+4VKsgyJXscOwbduunVN598oo/5+5u6zMw1ErOEEKKwJSfr4ebVq3ViU/ny8PXXZrcKIiKwBbfDkpLE3UDlss15befLVAmwcP1bdhWPUwQmXJ0LDgyUgh05UQq++gpef12/7lOmgB3FlQqDaT1hZw+D38zkdydEAf35J7RoAe+9pzdASMxmq8BClnjJyu5nxmOkJGEAyrBw25uPXQ3AmVkrBVDJGkdqsk2vFZaecPbGjdM1vps2hehoeOMNl9n0wpQg7O3tTXx8vASTfFBKER8fj3c+C8QLUawlJure0F136XoDixfDvHlww8oIM4SHw+P1d1P3n6Uow4Jyc8Pw9tJzl1mwVAvAgzRO7jgJ1arpJVQig9UK8fH653794Isv9KhHrVrmtusGpnwUCAgIIDY2ltOnT5tx+SLP29ubgPTNvIUQ9jMMWLpUZzq9/z6UKWN2i0g4mcx3vZYzaOUj1KzZiK0zI7mz0RVYu1bP92azIsSndiAAZ7bHEvDxxy6RROYydu+Gvn115nh4OFSuDP37m92qLJkShD08PKhRo4YZlxZCFDcJCbqy1Ntv6x7v1q1QsqTZrQJg7cSNVH6nHwPT9nClzy4Gfl6fEiWa6jvvuivH55ZpXosw2pF6yiIBOF1Kil5uNG4clC6tt5V08d+NZEcLIW5ev/0GDRrABx9krKN1gQB89qdVxPi15p632lDauMT+j5cy9Ov6eSp6det99bmfMHZ4tNDrmLp1g82bnddoV/fPP3D77fDOO3prx927oWdPCcJCCFHoTp3Sy40eeQT8/HQi1sMPm90qlILfh/+B75MdqHN2E8riToXfv6L2Kw/m+VylS0PZsleTog0Dfv1VB57iqmJF/Uv59Vf44Qd9uwiQICyEuPkMGACLFukCFlu36p2PTHYw8hz3t1dsnLweBTr7WdmIWRyW73MuS3uAR359Xq91heKXIR0eDp0764S7UqVg/XqX+LCVFxKEhRA3h6NH4cQJ/fMHH0BUlK4g5elparPSUhVL/jOf0i1qU3PjN+xoVZUUdw/SDAupbu6Mu1SRkKi4fJ27pFcqfmf3gbe3XudcXNYKnz8PL76oM8f37oUjR/RxFx96zoprLJQSQoj8stlg5kwYNgy6dNFDkTVrmt0qAHYtPUzCf/5L5wvL2HtLK853d2NXeV961hpPqyPRbKraiMhb63AkNIZuzfzzfP4r5QKp8s96fSMwsHj0hH/7Df77X/2B6/XX9VrvIryDlARhIUTRFROj14D+8Qfcf7/LbLiQHB7BvgFTuS1mMcpwI+q5T2g66yW2jFoOQKR/EJH+GbXfjyXkr1iItXIAlQ/EkXjZhk/Tproa1M3MZtOvsZ8fhIToRKwchETFMSU0hmMJiVTx9WFox7r5+rDjTBKEhRBF0+LFesMFHx9dbvLZZ3MdjiyMN+Wo6RHUG9Se+ioZw4DLc3+iWW+9T3kVXx/isgi4VXx98nUtt+qBeKxP4+j2k9T86qsCtdtlKaVHNzp00EPuCxfq77lMM4RExTFiYTSJqbo2eFxCIiMWRgO4VCCWOWEhRNGSXhmqVSu9BGXPHujTx64APGJhNHEJiSgy3pTzOx97o/Onkvm95bucHPguHioFN2xYLAal42KuPWZox7r4eLhlep6PhxtDO9bN1zXdb2/KN/TmWKytQG13WYcP68Srnj31Hs8AVarYNc8/JTTmWgBOl5hqZUpoTDbPMIcEYSFE0ZCYCMOHQ9u2OhCXLw9z5kClSnY93Zlvyuve38QJ/+Z03TaGqjXccfPxBDc3HSyCg689rlszfyZ2b4S/rw8G4O/rw8TujfLdMyvb6S768A1/J/nrTOGWLfWa4aLOZoPPPtNrvNeu1UU38rhNY3ZD/Pkd+ncWGY4WQri+dev03O/+/XqbweTkPBfgd8ab8sm/L7Gt00g6HfiUkx4B7Ju6mPqvdoaICB0Usyg72a2Zv8OGQ3X1WkXcwVSoBmzbBocOuVx95DwbPVpXverQQdd8rl49z6dw9NC/s0gQFkK4rkuXYOhQnf1csyasWgXt2+frVI54Uw6JimPprIXU2r2VLZZOHN9Sm9WXZ7Ct1UCaLplA5XKl9QNbt8625rMjXJvbPneFMzzL6u97Qa+X9Z1FNUM6JQXOnYNbb4WXXoI6deCZZ/K97Ghox7qZ5oShYEP/ziLD0UII1+XurgswvP663oIunwEYCj4fGxIVx/8+WcDHs9/g9XXz+Ca8L+4l4/n+i03cHvEpHukB2MkyzW0bBgluZeH4EX5L3w+nKAbhzZv1UPoTT+hErMqVoVevAq37dfTQv7NIT1gI4VpOntSVriZO1GUIt20DL68Cnzb9zTe/2dHvL4nhjaUr8LamYACeJPNgve+ZffZZ+hS4dfa7cW471rMyVZKO82b4YR6uUKFoFey4fFnXep42TQfeceMcWnDDkUP/ziJBWAjhGpTSe/u++qoehu7aFTp2dEgATpffN+U9oUd4e9wkHkxeiRUDZRikurnzZ9VGhZ7oc+P1jvtUpFVilD7eqRNUrVqo7cm3PXt05vPBg7rM6KRJuhh2MSNBWAhhvvQ34hUr9BZ+X34JQUG5P8/JkpJg4phkXpzciracZ1SD19ndpCKtYnfpalf+QfgXcqLPjXPbx0pVwP9sHJW9fODbbwu1LfmilO7tVq2q532//hruvdeUpqTPrcfFJ+Hv521KMQ8JwkII8738MmzcqJel/Pe/YDE/XWXz/w7S+53qxOzzwve+L6j6wq38svssialWIgMbAOYk+gztWJcRv+ygZtx+dleswXr/Flw+civP1LutUNuRZ0rBzz/Dp59CaKjeUnL5ctOaExIVx/Bfoqm26jQ9ow4S3fkWRlxOAQq3mIcEYSGEOaKjdfnBKlV08LVYXGIo9fKydRx54T2axq2lc/k5fBLamw4dugLg5gJlELuVuMS9i0Zwy+6/GPTwMNY2foKQiBdZXg7dqxw+HP7+W+8q5Cri4mDgQL3NYPPmcPq06a/1e98d4db5ikUnX8STFFJ/c6dnyXFMCfUs3NdUKeXUrxYtWighhLgmKUmpUaOUcndX6tlnzW5NJlEDZysrhlKgrIabuvLzYrOblCEpSakxY5Ty9FSqRAmlOnRQ6sgRdeigTZXjjJr7yXml5s1TCpTau9fs1mpWq1IzZypVpoxSPj5KTZmiVGqqqU1KTVVq8mSlDPc0NcLtPWXTfXSValjU5La9VfVhjnnNga3Kjhhp/piPEKL42LgRmjXTWbD/+Y/ecvCqkKg42kxaTY3hS2gzabXDykna49Qp+LHxBBp//gIGCtAdc599OwqtDbkaPRrGjIHHHtM93dBQCAykilc88ZSnwpK5eiclcJ0MaaV077xlSz3y8cYbeS6y4khRUXDHHTBqWAq31DnDzkfKkmpxx4pOtNtUtVGhF/OQICyEKBzffQd3362XpSxbppOIypcHnF/XOTtKwbffKIKC4H97GhMT9Ch4+2RZctIU8fE64IIOYMuXw/ff61KdW7ZAWBgelfxIwgsjLja9hJa5a4VTU2HKFD3k7OYGS5fqIiu3mTdnnV7x9Pbboerh9cRXacyy535lf/069Hh6Ih+27UXPHuPZU71hoc/xSxAWQjjXpUv6e6dOes/fXbv0z9cxo9j+kah4wvx7c6DPWOrVg/F/dSFo90KM1WF6nXJYmFOrXuVIKf0hpV496N1b3y5fXi/ZSjduHLz2GhgGp70C8D59FPyvzmWa1RPeskX3et98E378UR8rV86ha3/zau1aaNIEpk++wMo6LxFyti2lPZO5o3ktJnZvxMkGzZnR+klONmhuSjEPScwSQjjHyZPwyivwzz+6lrKfny7AkYXCLLZvXb+RA698QsXtodzLJUp3focxv12XkO3kkpO52r9fL9davVq3Y+bMrINYQICuJgacLxNI6Qux4O2tt3SsU6dw23x90Y1KlWDRIujWrXDbcIOEBP1ZIHp2BONLzeJh3yV47T0DQ4boD1mlStEN87c1lCAshHAspeCbb3Qv7fJlvfuNUjk+pbCK7R/6+FcCh3SnLjZsGJyd+g13vtrLodcokHXr9KYF3t4wYwb075/9cq3AQF1r+fJlkvwCqHBmvV6CO3duoTYZ0CMcn3/uMkU3Fi3Sydg1TkSw1r09HpeT9Vz/rFl6IxAXIsPRQgjHOXUKHngAnnsO6teH7dvh7bfBwyPHpzl6n90bJSfrjtqSV1diQe+9a7hZKJ/kInWWL17U3++4QweyPXv095zWS6fP/8bFcfje3ryn3ubcuav3paY6tbmAnq8+fFj/PGqUHvedMcPUAHz8ODz+OHTvrujr/g2L276Pp0rBUDb9uzx9OveTFDIJwkIIxylbVs8BT5+ue3V2Vr1yZrH9yO/2MC3wQ8aOhTMdel5LvDKclHiVpyzvc+fghRegUSP9e/P2zqijnJvrMqHV/Q/wFX31VPCwYXonImdRSieHBQXpbSVBD0G3beu8a9rRpDlz9Oe+nb8f5O/anRh7tA+3XInTCXaukmiXBRmOFkIUTFQUvPuurvtcurSe/81HIo6ji+1fjE9h7UOTeWDzOGpYStHyh96079EaIsKy3eu3oNKzvNOTzNKzvOGGuUeldOLSkCG6R/nqq3mvEta8OWzaBA0aUDUyiSbs5cSeajS55RYd3C9dcnzBjsOHdUWzZct0r33aNMeeP49CouJ477sj7PmpNilHbmGC/1SGJo/BctySUX3tzz+d9no7ggRhIUT+JCbq4PvBBzpzd98+aNHC1ExYACIiODRuHoSG0sX6D5G1n6Lu8o9pX7OCvt+JiVc5ZXlfC8IXLsCTT+p1vrffrr83bZr3i5UuDXfeCUD1pCi205zlq3+BttctU6pXryD/nMzWr4cHH9Q/T5sGgwbpHqZJftkSx39HXKRWuJVXmcnFBjaG7RrNibvbU+n7rzNGCsxOtMuFBGEhRN6tXq2Thv7+G/r21etCb7nF7FZxbmkEJR5uTzVrImBweOBkmn/2ZqFd364s79Kl9dDoJ5/ozesLEsh++glKlsSvxR0ApB6MhZ5N9H2OCsIpKbq9zZtDjx56jr9atYKf1w4h2ZQJ3b4dnu1WmubHDrCSDrgZaaTGuPPGg4OJuKcrG9IDcBEgQVgIkTdKwYQJuse7ejXcd5/ZLUIpWP7un/w5aTWjru73m2YYXDLO5frc62X3pm+v7LK82188rOdMf/hBr+X99VfHjBhMmgSVKuG2+EGSDO+rBTs66/sKulY4KUmvRV60CLZu1RsufPml3U8v6O8yq6H9Yf/bxbxPy/Drt6Vp7b6NXzy745mSgqEAaxoVL5/j2PmkvP5LTSWJWUKI3CkFCxboQvyGAfPnw44dLhGAj26PZ1VAHx58txUWz8ukurmTZlhIdXNn3KWKdlfdckTVrhuzvEslX2Fc2BfMnjFIr5dOzyZ21JB9QIDu8V4t2OFz5qg+NnhwwXrB6RUuxo/XUwwpKXl6uiN+lzcO7ScdKcc/s+5ixdcQWnsQa5Pbo9wNUq57vc0oO1lQ0hMWQuTsyBE9bLpkCQwdCu+/r7NhTWZNU6zot4AW37xCMGf5rPYA5ne5ne2nx9HqSLTe7/fWOhy5fj42B3bN5+Yi/XFTQmNosHk148K+oMLFeIyBA3VAK1Mm7//QnAQGXivYcaFMAGXOx4KXV/4Tpq5c0clis2dDjRp6vrpDhzyfxhG/y/QhfFuSO+fC63Hpr2p4lzlPVOnG1Io5yt//eZ4nAztT7djf117vPdUbMrGQy04WlARhIUTWrFa91Oitt8Bmg6lT9b6/LmDXLohp9xLdT81kX5mWpP24gg/XxqKASP8gIv0zlkbZW3XLUVW7rmV5P/sNVK8CsxZfS6ByuICAawU7NrR7hwU/GaywgSUlSZeMyuuHJS8vvUb5jTf0ZhElS+arWfb8LnMbrq7i68P+zWWos/wSfRPnsL5eI/Y9WJofj/Zg1OuPclurVrwdFceU0LLM8A+iiq8PE03YWrKgJAgLIbI2caJOwunYURdhqFHD7BaREr6B9e+tYcy69lQp8TDVnq5F87mDMTzcqfJXfIGqbhW4apfVqn9Pd9+ts50/+0yv+82lUEmBpCcgxcaSevd9hP2o66VUevEpPfwdHZ37OeLiYMQIneVesSKsWVPgnY5y+13mtpTrxAlQYa25bVUUYXTAk1RS97nzbKv3afjWK3A10Dp6WZsZZE5YCJEhKUkPP4NeY/n993pNqAsE4Jj3fsTtvrbct+YdVtGemT/60uK71zE8dMAoaNWtAj1/xw5o00aPFMybp4+VLu3cAAzQtaueE65dm9vKnKYzi4nbe1EH59x2UrLZ9EhHUJDOst68WR93wFaDuf0usxuufn95DF9/rYtuxK87ycKyvfEiFQNwUzZGlTpV5IPujSQICyG09GSc7t31G7Sfn97z1+R1vxfPprK49Xhqjn4GCzYsKLxI4Za/wjM9rqBVt/L1/CtX9B55LVronud332XaI9npypTR2dYWC7XjN7GYrlzYvFcPUyckZOxgdaPdu+Gee3SB5TvugJ07oUsXhzUrt99lVsPVqQk+RH3RmOefhzfLf8Ue94ZUTjmmP8i4ueHm7UWDpx92WBtdhQxHC1HcnTunt5v58kvd450wIe/Vm5zkj8+24/dab7qkRhPj34468Rt1XeRsShAWdHgyz8+fPh0mT9blG99/X39wKUzpy8WaNaNcEz00fXnvUWifMUydZZb0uHGwdy/Mnau3SnTCB62cfpfXD1crm8HFrdVJWF8Xi5vS+1aUKYFlflu9g1RcnEtXvCooCcJCFGc7d8L998OZMzrzefQNuC31AAAgAElEQVTofCfjONLp07qS467vFEvcz7NnUggxHVry4ayF1Nq9lQP1W/KQd1VM2Szv9Gk9ZN+iha4adeeduldpBsPQRT+6dcN33HgA0g7FQkAWBTvWrdM1pevWzcierljRhEbr4eoRC6M5H1eC+GWN4YQP48sNI/iJsrQe8Daop+A/T+l/X9WqN2XwTSdBWIjiyGrVlZpq19Zrfd98E5o1c+ol7SneoDZG8PdbX7I5wsYC9TVvjW6G3xsH2Lz/lE7kKVuTFa1rArA+q5rMzqSUnu997TVdHWzvXp14ZVYATnd1rbBR3o8kwxvLsaNQ/z+6itltt2VsrDt7tp5e+P5704Jvuk5B/sw/VoZjc3fwjmUoXUosJeDsMVD99QPMLn1aiCQIC1GcpC87mjlTF/8vXVpXccqFM6of3bixwal5ofg925laykpNDFrP7kmNfvcDHg5Zd1ogBw7orQXDwuCuu/S+tCbWTc4kIAAOHgTD4Ix3AD5nYqFCBXj9dfjlF50sdupUxrIjk61frzeO8o/5k/XGg1isaRiJBnz4of6AU8zkOvFjGEZZwzCWGYaxwjCMRYZheBqGccQwjPCrX40Ko6FCiAKKjtYZvK+8orNnL1+262nOqH4EGUHUmqZY1vcnfHo/jkXpxxhuFmqc3nLtsY5aw5sv0dF6q8EtW/QHmPXroUED51/XXtdlQn993zze8xirj69aped7K1fWbZ8yxdSphgsXdMJ927Z6f+e5nX/BTaVhgM5BSE42rW1msif7oicwVSnVATgBDAd+UEoFX/2yYyGaEMI0qam64Ebz5hkZvMuWZSrkkNMeuDkFUHtlFywPHXCjS+t42nzVl4QSVVCeXlnu9ZvdWl2nlig8e1Z/b9hQZ0Dv3q2jiIskrV0TEAAXL8KVK6Q0b8XGU7VIS0O/1pMm6aVHzZub2sTff9fLjn754gzzO81n504IHNlbD+e78F6/hSHXvyal1HSl1MqrNysAaUAXwzA2G4YxxzAMGdIWwpW5u+s34mee0dWQnn4605xbbj1dR/RC/xUs0xRtfj3M8a/vZsvB8qwft46AhF1YwtfA2LF62Pe6ZJyCrgHOk4sX9WjBbbdl1MoePVovBXJFQ4bo9d0lStDIfQ/9bTM4HmuFF1/U/w4HrPvNr5Mn4amn4OGHFc8Y3xFXNoiea/pR8sJx/fquXp3l612c2P3qGIbRGrgFWAl8rZQ6bhjGt8BDwG83PLY/0B+gatWqjmutEMI+8fEwcqT+CgzUdZ+9vLJ8aG7zrQWuJEVGNmzQoZ102LaVO/ftoZl1B4HNK/LassepWPHqfrrZ7P16fU3m/M5L22XxYl0nOzZWfy9d2rHndwZv72s/Nji1hicZyJYd3QisXtm0JikF334L378cwQOXF/FhwHoCYjdBq1Z6KVzlq21zhb1+Q0J0wZK65tSctisIG4ZRDvgUeAw4oZRKH7zfCtS+8fFKqVnALICWLVsqxzRVCJErpXT265AhOiu2bVvd880mAEPuPd30AHp9oM5rL7RbM39KRvxJ8PgRuKs0ANbeP5hJKx4DOxNhnVqiMC0NevbUO0U1aAAbNpgfHOx1/jwMGwaPP06pegEAJEQfhYfNCcIHD+pO+MWVEayxtMfLlogRi/6b/OAD10loA12UZsIEXWZ01ixTmmBPYpYn8BMwQil1GJhnGEYTwzDcgG7AX05uoxDCHv/8A5066WHn226DyEgdgHOR23xrQStRgd6Mp/rgkXikJ+K4uXFvu1tdZymKuzv4+uoiFpGRRScAg/6A9cUXEBFxrWDHlX25lKx0AqsVPvpIT6HHbjjM1K5r8DJSrr3eVKzoGgHYZtO98TNn9Px+SAh8/rlpzbEnw6Av0BwYaRhGOLALmAdsByKUUquc1zwhhN0mT4aNG+HTT3VPrpF9CxfsmW/t1syfDcPbcXBSZzYMb2d3AI4/eoXneybTqRP8VG4ANg/PLBOvTBETowuVREXp2198oYfvPT3NbVdeeXvrJUmxsZQO0j3htEOFG4Sjo/XKreGvJTPL/112pdam9Z1Kv86ulHi1Zw/ce69eIzVnjj5WpYrza3znRCnl1K8WLVooIYSTbN6s1F9/6Z/j45U6ciRfp1kUGavumhimqg9brO6aGKYWRcYWqFk2m1JhI8PUQUtN9a4xWr39tlJJSUqpjRuVmjBBfzdLcrJS772nlKenUr6+Sv36q3ltcZRmzZR66CGlbDaVZHipRbXeKJTLJiUpNWqUUu7uSnUqu1El+NdXCpT6z3+UOnnSNV7v9IaOHq2Uh4dS5cop9fXX+o/UiYCtyo4YKUFYiKLowgWlBg9WymJRqksXs1uTydEdZ1VoYF+lQB3xqqX+nrPG7CZl2LhRqQYN9FvfU08pdfy42S1yjIcfVqpRI6WUUn3b7FH3Njvv9EtO+OqUutc3TA1ngppX7j/KZhhKBQYqtXix06+dZ6+8ol/zp5/WHw4Kgb1BWJYXCVHU/P673v0mNlavW50wwewWAXqqbcNTn9D055G04wpb2g2j+a+jcSvlxLW8ebV8ua4a8fvvDt01yHQ1asChQwAYQfXY+7vzLnXhAjzV7xIJPx0gjC54kowtwUJo3btImz2HLnebk2X8L+fP6+VmAQG6bGenTvDgg2a36l9cbNW5ECJHCxbAww9D2bJ63vfzz/XPJtu7F3rWW0vLn9+kFJewurlzpdddrhGAf/tNV48CXbRk166bKwCD3pDhL50j2zY1jD4nJzmlANWSJTp5fMtPicz06Yc3Sbhjw6Js7Kh4GxP/sL+KmlMtWqSrg/Tpo2/7+7tkAAYJwkK4PpvtWi+HRx6Bzz6DbdtcIoM3JVmxpNtsVjQYTI2D6/FI34Ddlsa2bxflqbSlwx0/Dk88oX9nH32kj3l5FY21vwXQ5EwYY3mbuCPW3B9sp1On9N4PXboo/qO+J8arDg2S9mK1WEgzLKS6ubOpaqPCKSOak2PH9H7Y3bvrZLWJE81tjx0kCAvhynbuhLvv1jv1XLqkg8jAgS6RwfvXz/vZ7teOzr/2p5nPNvZ19SbV3f3am/If/g3yVNrSYWw2vWNQUJAedh43TveMbmb790PnzrBpE541A/AgjZPRpwp82vSNo4KC4M+fj7Kvdhfej+vJiYqV6fzcJzz19CSm3vMMPXuMJ9I/yLllRHOzYYNu6LJlulznli1w++3mtcdOMicshCtKTNTB4/339XDz1Kkusc8vwOWEVNZ0nUr7P8aQYnjx16DZ9C5RCWUY9OwxnlZHotlUtRGR/kEYZvSMliyB/v31UpRZs6BOncJvQ2GzWGDpUnjiCUrX12uFE6KPQvf8F+w4fFgX3QgN1YMu3w2KpcaA9TBtGvvbPMqRX3eTmGol0j8IcGIZ0dykpel13k2a6Kma0aOhVq3Cb0c+SRAWwtUcP64rXR04AM8+q6sMlS9vdqsgIoL9s8OZtLgBH50eT0yNB7lt2Wc0qVuFKpNWE5eQSKR/0LU3ZXDyBgvXS02FHTugRQs93/vrr/q7q2224Czpda2PHqVc+yZAesGOO/J8KqtVz3j8OjyC7mkLmBhs0HjVVNzcWsNDR8DXl24Abm7OLyOak5QUvTb+55/hzz+hVCndbS9iJAgL4SrSP9FXqqQLG8ycCe3bm90qAM7/Gk7J7h2oYbPxmeHJwXHzaTqy+7X7HVHaMt82b4Z+/fS8+cGD4Oene0TFyXUFO3zq6J5w2uG8z8fv2gV9+4LHn+tYY7TXWw2GA6s6QseOuqrYVU4tI5qbTZv0a75rF/TooTewuK6GdlFSTD4mCuHClIJvvoHateHoUV3KcfZslwjASsGad9ZgffQx3G2puGPF25JCfUvmuV5HlLbMs0uX4NVX9VhpfLzuBfn5Oe96ri4wUP/9+PlxT8NzzCs7yO6nJifrUdxmzaDs3j9ZUeZJ3K8rMUpkpNOanSeJiXpnqLvu0mulFi+GH37I9OGgqJGesBBm2rcPBgyANWv0G4sLbWwetyuBnQ8NpeORLznp4Y8NLyy2tGxLThZqz+j8eV10/9AhvVZ64kSXWKplqhYt9N+PYXBLDV8OH7bvaRs36k7lnj3Q/8kEZiy5H4uPt07+s1pdp+Qk6LZs3gyDBsH48TdHprs9FT0K8iUVs4TIgs2m1NixSnl5KVW2rFIzZypltZrdKqWUbsb0z21qs+VOlYqb2nzfmyr1/GXXKEGYlJTx85gxSq1fb15bXNi84C/VeO/3cnzMhQtKDRqklGEo1fnWLWrpkqtlHFeuVOr8edd4vZVS6sQJpfr3V+r0aX07Odnc9tgJKVsphIt7/nnXKp24caM6/cIINbBRuAKlXmsRro6GbDW7VZrNptS8eUpVqqRUZKTZrXF50c16qUNUVZcuZX3/kiVKPVJxo3qPUWpXtQd1KAgJKdxG5sZmU2ruXF3r2dPT9dqXCwnCQrias2eVeuEFpbZeDWypqea25zqpazeoNIuHsoFKwV0tHrnR2fXt7XfokFKdOum3q1atlNq92+wWuabVq5Vq3FipfftUdNcRKgV3tWdnWqaHnDqlyye3YoNKNjyVTU/7K9Wnj1KJiSY1PAt//63U/ffrtrVpUyRfc3uDsCRmCeFsSunkkXr14KuvdGYn6ExoF7Bj4QGO398LN5uuduVuUXQuGe4aW/1On67rJK5fD598An/8oQsyiH+zWPQyrcOH8bpNF+w4FX0S0H+C8+frX91PP8HPlQfjqa7b67dOHdfKLh41Si87mj4d1q27qV9zCcJCONPff+vC8U8/DdWqwdatuuKVC7h8GeY/9D21H2tEudQT2Nw99F6/Xi6UiHPqlF4zvXs3vPyya2wK76oC9F7CHD3K36VLATD2k9W0GLaRlvck0buXjbq1rERFgf/L3fUeuq601+9ff+m18aCL0+zerZPubvK13q7xUVyIm9UPP0BEhO7FvfSSywSRVSsV/V80KHOwPvVrdKbW0k+wnDsM4eH6DdmsutTJybpSWKtWugzj22/rN2GX6Ja7uKsFO/Zs2c0XSdW5m5J4xXgTtfUO6qq97PV/gVqPdsfSYCg0GKFfZ7Nfb9DLjsaO1dXhunXTxTcqVTKvPYVMgrAQjrZxo35jad8ehg6F557LqGhksrNxiUR0epcTO0/jUWcOn6xtSvO2P1+9t4q5b8YbNui1Mnv36t9b584u86GlSLhasCNm2x6igu+lk/digs/+QbdSz9Hryv9IPuuNJeC6tcOtW5u/CcjatfDCC7r29XPP6epwxczN3c8XojCdO6fX/LZpoysfKKU3XHCBAKwUbOs3HRUYSOedk2nSGP6KtNK2rdktQxddGDhQb1SRmKgL8L//vtmtKpq6dmWvTwWaH9vLqpSOjOVt+l2az5bA+rTvOx169jS7hRl++kn3wq1WvdXkV19BuXJmt6rQSRAWoqCUgh9/1Mkjs2fDa6/pzeNdZAj12J7zbK7yCC3mDKScisfm4UXLmf3wLukivczffoMZM2DwYL1rVKdOZreo6Jozh9879aLVkWg8bTrxyorBxmpN8PSvYnbrtHPn9PcHH4R334XoaJeoDmcWCcJCFNScOXqz1cBAnXj14Ye6mLzJbDZdfrrdHZdoeGIlCgMDsNjS9FygmU6dghUr9M89e+qs3mnTXOL3VtQN7ViXqJpNSXb3JM2wkOLuQVTNpubscHS9Eyf0/s6tW+taz6VKwTvvQIkS5rbLZDInLERB9ekDzZvrrdRcZA7zwB8nWPfMLP57+G3atfPnfK8QSr7UTe88Y2Y2bPoGta++qkcKjhzRb8ING5rTnpvNnDl0e+MN3Bb/yWBvd2rt3sqB+i15qn938zZbUAq+/hpef11PN7zzjsv8P3EFEoSFKCh3dx2EXUBqiiK0x9e0WfQ6T5NIuTGP8Mg7TTCMDlA3zNxs2EOH9Jx5aKiukz17drHvBTmctzckJNDVz0bXGS+b3Ro4exaefBLCwuCee/RrXtfkHrmLkSAsREFcuAAjRugMz6ZNzWtHRATHpy8k/pdwuiRuZXf5tlQImU23NtdtaG9mNuzx47q3axjw6ad6udZNvv7TFOlrhWNjdXEYs5Utq1/nGTOgf395zbMgQViIgti+XVf16dzZtCCctCYCS4f2VEpLpBKwr9ub1P9lYqY3vJCoOHM2YI+P19sLVq4MkybpfX6rVnX+dYurQL2XMEeP5vgwp/497NgBw4bBt9/qPY5DQ10mSdEVyccSIQoiKkp/N2k4etPsaKZ3X4UlLaMEYZ07fP8VgEcsjCYuIREFxCUkMmJhNCFRed/03W4pKTrztWrVjL1oBw2SAOxs6cvhYmOzfYjT/h6SknS5yRYtYNs2vfYXJADnQoKwEAURGamr+xRyhZ9zx5NY0uQtWvZvRmVbrC416eaW5V6/U0JjSEy1ZjqWmGplSmiMcxr355/6Q8mYMboCUnrvTDifl5cueJJDrWWn/D388Qc0a6b3+H36ab058V135f98xYgMRwtRAOc3bmZ3mao8PXxJoQ3zrh27Dv93X6CzdR9bG/Wh27KJuB3pk23S1bGExCzPk93xAhk+XBfa8PeHxYv1ML0oXLNn53i3U/4ePv5YZz4vXw4dO+b/PMWQ9ISFyKeQrUc4fyaBrX41nDrMGxIVR///fsqkVn2ZVa4f975zL95uqRyYsZKWO74m9FQibdYmUuN8Y9qsTfzX9av4+mR53uyOF0jJkrro/q5dEoDNlJSU7V0O+3tYtgz27dM/z5ihC61IAM4zCcJC5NOUVQdo++KXfHT309eOOXqYNyQqjh8/XsDHs4byxp/f0OvcfH6p2JmtK0OpNeB+u+b3hnasi49H5nWZPh5ujinecPasrvm7eLG+PWoUfP45lClT8HOL/HnjDZ0Il40C/z2cOQO9esFDD8HkyfpY+fJSaCWfJAgLkU/pw3c2i1uWxx3hs1nbGfn9bLxsybhjxd1I5e96fkzZqLNf7Znf69bMn4ndG+Hv64MB+Pv6MLF7o4IPm//yC9Svr4tvxFy9niThmK9CBUhIgEuXsrzbnr+HkKg42kxaTY3hS2gzabX+UKcU/O9/+jX/8UdddGP69EL6R928ZE5YiHwaufUnypyK482HhmQ67ohh3rRUxYqn57Lg59cpxUWshjs2bKS5ubOpaqNrgd7e+b1uzfwdN1d9/LjOdF64UCfjLFumvwvXYMda4Zz+HtJHV9I/3KWPrlRb9D3Nxr4Jt9+ui280auSU5hc3EoSFyKfuJ3dw4EJapmOOGObduewol556nocuriKixJ0M79YfP8tZWh2JZlPVRkT6B+F/NdBX8fUhLotA7JT53nQrV8KSJXrd7+uv64phwnVcv1Y4HwU7rh9dMZSNCpfOcaq0H2+UqE/Y9Om66IaUnXQYGY4WIj/S0ij3917K3dPKYcO8iYm6xsFDXSwEXo4h6oXpnFz3E6dqViHSP4jprZ8k0j8oU6B36nzv9Q4ehN9/1z/36qWHn4cNkwDsiq7vCedD+ihK9bNx/PDDW/z4wwi80lL45wo66U4CsEPJ/yAh8iMmBhITqdXxHjb0alfg0+1+ax5JU2ewPnkKHfq2ocT4AzS71ZNmABZLttWN0r87rfqR1QqffQZvvQW+vtChg16LWq2aY84vHM/fXydn1a+fr6cHlvak06ofeO2P70hx82DcfX1JdvO4NvoiHMtQSjn1Ai1btlRbt2516jWEKHTz5+seYXR0gXYASjiRxK57XuSuA98CYPP0xi18tXk1nq+3ezf07QubNulM2JkzpfDGze7kSc6168gtu/8itHYr3n7gv5wq7YePh5tjkvmKEcMwtimlWub2OOkJC5Ef3t7Qpk2BiuSvG7+eKqNfoI01BgUYgJs1VRfdMDsIHz2qk61Kl9YfOJ5+WjKfi5LLl3WGtH8eg2b58txSqxqbew/gPVttTp9Pwr8wa40XQ9ITFqKQHT8Og19KZVJIXbw8FdZBg6k6462MvX7Dwszd7Sh9jemXX+oNFypWNKctIv+6dYN//tGbKeQmIkJXOvv5Z728STiEvT1hScwSIq+U0nOl+Xja8tdW0Cwoid+WebB6yO9UPLmTqh8O0YF37FjzAnB6Vlj16rr4PugaxBKAi6aAgFx3UuLyZRgyRI/oHDwIhw8XTttEJhKEhcirQ4d0RaiFC+1+yrEvl/BP6cZ0+qgjo/xmsGMH9PuoAR6+JXVhhBzKTjrdunXQpImu+dy7N9x2W+FeXzheYGCOBTtYtUrnMnz8sc543rkTWubaaRNOIHPCQuRVVBRcuZKxFCQHaamKzZ3eofXqcQBYLe68NKcFljr6/uwKIwCFMwf36qswbRrUrKl74e0KnuktXED632ZcHNTNYrnazJng4aE/gN1zT+G2TWQiPWEh8ioqSq+VzKVi0Pbt8GPAG9x1NQAbgJuhsERsuPaYQt9m8Ebly+tAvGOHBOCbyfUFO9L9+mtGedFZs+CvvyQAuwAJwkLkVVSUzor2yXrdZOIlK6Nfv0TLljAn7VliurymH+vmphOvrtvvt1C3GQSIj9dDzr/9pm+PHAlTp+rdj8TNo359PdRcqxacOgU9euhkrQ8+0PeXK5ft368oXBKEhcirqKhsayVv/TqafeXv4rapL/Hss/DL/sbU/f1DjGwSrwp1m8H0DRd++EFnzoqbV/ny8PLLsGGDfs0XLdJ/f59/bnbLxA1kTliIvLBade3c64NwRARJi1ex7ed/uGPffC5YfHF7awi9x1/3vNats8x6HtqxbqY5YXBC2ckTJ2DgQJ1I1qIFrFihE7HEzW3ZMl1spVkzmDMn3xW0hHNJEBYiL9zcYPTojNsREdiC2+GVkkQbYG/1TlRdO4+GVcvbdTqnl50E2XChuDp2TE81vPii1Ht2YfK/UYi8OHpUb15+yy36dng4pKRgAMpioV7/tnBDAA6JissxyDp0m8F0R47okpqdO8Mzz0DbtlLvubjp18/sFgg7yJywEHnx+uuZ11MGB5Ni8SINNwwvr0xJV5CxBCkuIRFFxhIkp60Fttlgxgxo0EC/CScl6XKTEoCFcEkShIXIixuSsmy338kz7v9jWeusq10V6hKkAwf0MqOXXoJWrXQ5Qm9vx19HCOEwEoSFsNf58zrQXReEj645wM8pD1OmfkCWiVeFtgTp2DGdbLV9u07CWbFCl6AUQrg0CcJC2Ouvv/T35s2vHTq+NAqACvdlvZ2h05cgxcdfPWEVnYSzezc8/7zseCREESFBWAh7RemAe31POOXPKFLwoGbXBlk+ZWjHuvh4ZM5MdcgSpNRUGD8eqlaF9F3KXnxRB2MhRJEh2dFC2KtzZyhbFipVunao1IEoDng3pH4Zzyyf4pQlSFFRure7fTs89ZQOxEKIIkmCsBD2qlVLf6VTimrxkUTXeJicyiA4dAnSe+/prwoVdPGNRx91zHmFEKaQ4Wgh7JGUBN9/r6tPXRV/2sYQ21SOd+hTeO2wWPS63127JAALcROQICyEPaKjoWdPXYv3qr92ujGfXlTo7sSdaK5cgTfegN9/17dHjoS5c3UBfiFEkSdBWAh7ZJGUFbdkO03Y7rwyzOvW6WVHH34ImzfrY5L1LMRNRYKwEPaIitJJWTVqXDsU9NO7/OzegwoVHHytS5dg0CC49169YUT6DkxCiJtOrkHYMIyyhmEsMwxjhWEYiwzD8DQMY45hGBGGYYwqjEYKYbqoKGjaNFNPtMqJSOIqZL2lYYEsXgzTp8PgwXoYvF07x19DCOES7OkJ9wSmKqU6ACeAHoCbUqo1UNMwjNrObKAQprNaYceOTEPRycfiqZJ6hCv1mufwxDw4fx5Wr9Y/P/WULgwybRqULOmY8wshXFKuS5SUUtOvu1kBeAaYdvX2CuBuYL/jmyaEi7BYICYmUy/46O/bqQX43OWAnvDSpXqP4osX9S5NZcpAo0YFP68QwuXZPSdsGEZr4BbgKJC+BcxZ4NYsHtvfMIythmFsPX36tEMaKoRpDAMCAyEg4Nqhc6sjAQjoWoAgfPYsPPusLgLi6wurVukALIQoNuwKwoZhlAM+BZ4HLgHphW9LZXUOpdQspVRLpVTLCg7PWhGikP3wA3z+eaZDi27pSyevNdRo6Ze/cyYkQMOG8N13MGoUbNsGt9/ugMYKIYoSexKzPIGfgBFKqcPANvQQNEAT4JDTWieE2bZvh7fegm++yXQ4IqYc55sF4+aWzfOyk5Skv/v66sSrLVt05rOXl2PaK4QoUuzpCfcFmgMjDcMIBwygl2EYU4EngSXOa54QJklJgdGjde80MREmT752l7p0mfabxtOh+j77z6cULFiglzht2aKPDRuWKdlLCFH85BqElVIzlFK3KKWCr359AwQDm4D7lFLnnd1IIQpd5866RnOPHrpE5H33Xbvr5Iq/GJU0ila+e+0718mT8PjjOus5IEAynoUQ1+SrWIdS6pxSaoFS6kTujxaiiEhO1suRAF55BX77DebNA7/M876nQnX1rFs72dGL/fFHqF8fliyBSZMgIkLfFkIIpGKWENrmzdC8OXz0kb7dtav+ykLa1ijO4Efd9gFZ3p/J339DnTp6bnnYMHCXjcuEEBkkCIviLTER3nwTWreGCxfsWp9b9p8o9pZoTslSWdRxVkoncS1dqm8PGwZ//AH16jm44UKIm4EEYVF8bd6sS1FOmQJ9+8LOndCxY87PsVqpcOFv4gOzGIqOjYUuXaBPH73TEeieb55TqIUQxYWMjYniKy1Nf61cCfffb9dTzl9yw892hklPJWYcVAq++gpee02f7+OP9QYMQgiRCwnCongJD4dNm2D4cLjrLti7Fzw87H76jh1gxZ36d5bOODh7Nrz4IgQHw5dfwm23ObzZQoibkwxHi+Lh4kUYOFAvNfrqK7h8WR/PQwAGUB9/wiSG0bTpdQf79dND2WFhEoCFEHkiQVjc/Fau1AlXM2bAq6/qTOV8rtWt9MdPBL4YRbMAACAASURBVLtvoHLl6w5aLNCggf4uhBB5IO8a4uZ25gx06wbe3jpLeepUKFEif+ey2fA/vZ3jlZplbKh0/rzeASky0mFNFkIUHxKExc1p82adMFW+PCxfDlFReg64ANL2HqCk7RIpDa7LjN6+Xc8JnzxZwAYLIYojCcLi5hIfD716wZ136opXAPfcAz4+OT/PDseX6kpZJe++LghH6WM0b17g8wshih/JjhY3j4UL4aWXdCB+5x3o1Mmhpz+6P4kUalK9c4OMg5GRULky3PqvbbWFECJX0hMWN4eBA+Gxx6BKFdi6Fd591+HbA4aUfZYGXn9Tp6FnxsGoKNkJSQiRb9ITFkWXUvrLYoH27cHfH4YOzfOyI3tt366ToK+d3mrVWx62aOGU6wkhbn4ShEXRdOwY/Pe/0KaNrv3cvbtTL6eOxjI7rC2L2n8GPKQPurlBTIz+ICCEEPkgw9GiaEkvEVm/vl7/m9/lRnl0bnUU1WwHCWhQ9t93Glls5CCEEHaQnrAoOg4fhhde0MH33nt1ichatbJ9eEhUHFNCYziWkEgVXx+GdqxLt2b++bp0/KoofDGo1LFJxsG334ZDh/Sew0IIkQ8ShEXRceyYXv87fbqu1ZxDhaqQqDhGLIwmMdUKQFxCIiMWRgPkKxCryEj2UYdGrUtlHFyxotB64kKIm5MMRwvXtn8/fP65/rl1azhyRM8F51IickpozLUAnC4x1cqU0Jh8NeOWQ1HsK9mcsumj0WlpejcHWR8shCgA6QkL15SWBh99pNf7+vhAjx7g5wdlytj19GMJiXYdt2vIOi2NMI8HOV7rnoxje/dCUpIsTxJCFIj0hIXriY7Wvd4334QOHfQORX5+eTpFFd+sK2Rdfzx9yDouIRFFxpB1SFRcpudcTnbn6QszudC1Z8bB9EpZEoSFEAUgQVi4losXdZnJw4fhxx8hJEQX4MijoR3r4uPhlumYj4cbQzvWvXbb3iHrI1+tZLgaz33eERkHS5XS2yLWrYsQQuSXDEcL17BvH9SuDaVLww8/wO23680X8il9SDmnoWa7hqw3bqTu4AcZjxU1xgfahule+qOP6i8hhCgACcLCXImJMGYMfPCBDr5PPgkPPuiQU3dr5p9jJnQVXx/isgjE14asY2Ohb18s6mpvOTUFwsOhVSs9Z+2kylxCiOJDhqOFedavh6ZN4f334bnn9PxvIcp2yLpDHZg5E1W/PmkH/iEFD9Jww/D0hOBgOHhQ99h/+aVQ2yuEuPlIEBbmGD0a2rbVtZdXrtSFN3x9C7UJ3Zr5M7F7I/x9fTAAf18fJnZvRLfmAVxevp7tnndQJ20PL9Vfy7nXxkLY1aHoqChIToaqVQu1vUKIm48MR4vCpZQu89ikCbzyCowfr5OcsuDIilfZuTZknZoKH3yA1VKRTz/1572Vs0nEhwkfGwwcWBM3t9YZT4qK0nWjGzVyaFuEEMWPBGFROM6ehdde0zWf0zdcyGHTBUdXvMrRV1/ByJFw4gRzpyfzSmxDOnUqwcyZUK1aFo+PioKgIPD2dmw7hBDFjgxHC+dbuFAH3/nz9TCuHRxd8SpLV67A00+j+vZFnThBMp78dKEj8+bB0qXZBGCQPYSFEA4jPWHhPCdPwqBB8PPPOmgtX64Tsexgb8WrAvnkE/jhBxQGFhTuhpWfB4VT6pnW2T/HatX/piZNsn+MEELYSYKwcJ6DB2HJEpgwAd54I09LenJdPkQ+54zPnoXYWC7VbMyYo4O5yC1MM17Fy0jBzcuTUl2Cc36+mxu89Zbd/w4hhMiJBGHhWIcPw7JlMGCAXk975Ei+im4M7Vg305wwZK54lec5Y6V0j/zll7nsXoZGlj0cjvXhpYEvYuvWGMuWcL38qHXmXvCNgX50Qx863HEbVKyY53+TEELcSIKwcAybTW8xOHy43uGoe3cdqPJZ9Sq3ilc5zRlnCsIREfDbb7BhA6xfz8FyzXn05By867mxfj20aQPQGu7/9xB0VoHe9uZoLl88Ssmjh/L17xJCiOtJEBYFt3cv9OunA13HjvDFFw7pKeZU8cquOeOICGjXDpWUBMCXXgN5+fw03nzbnZEjwcsr5+tnFeiDjh/gz4A6tLP/nyGEENmSICwK5vJl3Z1UCr75Bnr10uuAnSzXOePERAgPR6WkYgBpuJFawZ8tS93tXt57Y6AvnXyZagknWNC4gwRhIYRDyBIlkT/79+vAW7IkfPst7NkDvXsXSgCG7EtOvtmuJkyahKpZk98PNSTR5kkqbigPT178IThP9TVu3A6x/sl/ADheM6jA7RdCCJAgLPIqMVHP+wYFwYIF+ljnznDrrYXajKxKTk6vp3jkxUdhxAjCU9rQb9YdvHVHGJeGjsVjbRhud+ew9CgLNwb6BleD8AM9OznynyKEKMZkOFrYb906Pfe7fz88/3yhb7hwo27N/OmWdATWbILoXaiRP3Kx5K08776Q1TzKR99Ar163Yhh5C77Xnx8yksN2tLiXre0a8eD99q11FkKI3EgQFvZ55x0YOxZq1NAbLtx/v9kt0olX7dtDSgpKwZqSXeh+8Rs6PeXL7o8d0znPbTtEIYT4f3v3HR9Vlf5x/HMyJCTSQSyEpquyFAsKAqKSAIIoYhZFbOjaUATEAkgUERcRXYplVRTRFVmKIggWFJQQUYxSDFXwJyKowYIoTQhp5/fHSWimTMrMnUm+79crryTDnZnnziXzzL3nnOcpDV2OlsJZ676fdRbccw+sXRsaCfiPP2DQIOyBA5CdTVYOfGHaMvXtmsycGYCr4+np8OqrkJZWxg8sIhWZkrDkb/t2uO46eOIJ93vPnjBhgpuI5bXZs6FZM3JWrCQzx+cmXlWK4q7ZcVx2WYCec80a1/P4iy8C9AQiUhEpCcuRrIVp09zEq1mzDp0Jh4Jt29yHgSuvZEvGibSyy7mh4cf82HcUUUsWUSWfghtlJjXVfVfjBhEpQxoTlkN++MGVm5w/H9q0gZdfhubNvY7qIPvCi2S/+z6PVXmCx3bdy+AHKzF8OERHBzD55klNhZo1oXHjwD+XiFQYSsJyyLZt8Mkn8NRTrlOQz1f0fUrBrwYMb7wBixfz68U30G/FMFZn9qH2WaewbDKccUZAwztSaqobFw/SOmgRqRiUhCu6jRthwQIYNMid/f7wA9SoEfCnLbIBQ2amm3g1cSIAVV+Ywo6oRdw5vh2DBgX888GRsrLcmHC/fkF8UhGpCJSEK6rMTPj3v+Ff/4Jq1eD666FOnaAkYCiiAYP9xa1Hzh2HNUAUGcy5K5na9wbh0vPRKlWCLVtcL2ERkTKkiVkV0fLlcM45MHw4JCTA+vUuAQdRQQ0YjvtqFfbcc9nzzU886HuM/cSQE+HDFxNF7Z5xQY3xCMcfD/Xqeff8IlIuKQlXNLt2uXW+O3bA3Lnw+utBLzkJf63LXHvfLgA2V23PU3Ufo+Her9h0RSLp7y4i4tFRmEWL/tLrN2imTnXj5CIiZczYAC9BadWqlV2xYkVAn0P8sHIlnH22m1iUlOTOhIN06Tk/eWPC5/zfcoYseY2/7Ujj3BaL+Xr1uZxwgmHiROjRw7PwDvnsM7j6ancW/PnnXkcjImHCGLPSWtuqqO10Jlze/f67KzLRqhXMmeNu69jR0wQMbvLVm38kM/WNEZzx8yYqZ2ZyXOo++vY1fPVVCCTgP/6A2293bRqNgXHjPA5IRMojTcwqr6x1xTYGDnSXnhMTXbejUJCZCVdfTfM5c7C4iVcGy39v+5yTJsZ7HZ177Tp2dDOi77sPRo6EqlW9jkpEyiEl4fLqjjtg0iR32XnhQjjzTK8jOshWiuS7nbX5OPoOeqdPobLJwBcdxUk3xXkb2HffQf36EBkJY8e6yWqqkCUiAaTL0eVJTo5b0wpwySUukXz+eWgk4E2boGtXfv5wLZf3sPwtaRLPNZ9I2pRF+EZ7PPEqMxMefxyaNTs0AatzZyVgEQk4nQmXFxs3wm23ueSbmAiXX+51RK7VYFISbNuGfeUVMoji7iWb+YgWjBtnGDQIKlVqB3iUfPNi7NsX1q1zdamvuca7WESkwlESDncZGa7oxqhRrsPR7bd7HZGTkgLx8XDgAACrq1/Apbtn0LRTLGtfhL/9zeP4wL1uw4a5S9Dz5oXAbDARqWiUhMNZairceKPr8du7NzxdRp3sS8taWLwYm5GBAbKJ4O0DF/PoK7H8858el1+21l1+joqCCy905TpHjdLEKxHxhJJwOMvMhN274e23CVwj3WJKTgZgwwnxnEQ0PjLI8UXR/8146nT3NjS2bIH+/aFhQ5g4Edq2dV8iIh7xKwkbY44H3rTWXmCMiQW+ADbl/nMva+32QAUoR3n/fddYfuRIOPdc+OYbN5vXazt3wtCh8NJLfNOwIy1+/IjzYt6na+M3+P7Cv3NJbEMSvIotK8tdJRgxwp2GP/qoV5GIiByhyCRsjKkFTAGq5N7UBhhtrZ0YyMDkKNu3wz33wLRpbhbv0KFwzDGhkYDnzoU778T+8guTawxm0PePUL3lj2y9MJ3J0ZcA8MnhHZKCaf166NPHXbq/7DJ49ll3JiwiEgL8WaKUDfQGduf+3ha41RjzpTHmsYBFJo61rnZx06aut+7DD8OXX7oEHEBzU9No/3gSJw17j/aPJzE3Ne2vG6WkwM03wz/+wQ/pdWmd8wVjjxvLybdtpEaXNUREZx3cNK9DUtAdcwzs2QNvvukmXykBi0gIKTIJW2t3W2t3HXbT+0Ac0BpoZ4wJZmv1iufnn13hjdNOc2dzI0dC5coBfcq8us5pO/djOdTr92AithZefx3bqRM5U14jgyiu3/ksFw1rxerV8Gftn/J93II6J5W5efNcK0Rr4aST3PKtK67weEaYiMhflaRYx2fW2j3W2mwgFTj16A2MMX2NMSuMMSu2b9dwcbFlZbmSk9bCiSe6JgKffgrNmwfl6Qvr9cu330LnzuTccCM5+zOIyMkmgmxm9P+UMWMgJuavHZLyFHR7mUlLc2t9ExJg2TJXNxvA5wvs84qIlFBJkvACY8yJxphjgC7AuqM3sNZOsta2sta2qlu3bqmDrFBWrXIzdq+6yhW6AFfxKiJ4xc3yO2P15WRz6YJp2NNPJ+OzFTzDIA4QRbZxvX7rXRt3cNshXZsQE3lk4ouJ9DGka5PABJydDc895y7Zv/8+jBnjukYFuUeyiEhxlWSJ0iPAYiADeMFa68FAXzm0bx888giMHw/HHuv6/Hbs6Eko9WrGkHZYIo7OTOf16Ymc+fM3fFqnB713PE+T+Fiu7JdA/U3JEBd3RMnJvMlXYxd8zbad+6lXM4YhXZsEblLW3r1uxnPbtm7pUUhUAhERKZr6CYeKDh1gyRK45RZX87lWLc9CyRsTbvrdWtr+sI6U+mfQdclalvzUjflVrmL8BMPNN3s8xLp/v0u4Awe6GeI//OAqX2ncV0RCgL/9hFWsw0s7dkD16i6JPPCAm3QV730rv4SWsTSe9RpnTHsQCxwgmk4sol7Pdmx41g1Te+rDD91ktc2boUkT16KxQQOPgxIRKT51UfKCtTBjhhvDHDvW3da1a0gkYHbtgn79OGvMAxgsPiyRZPDSdcnMnu1xAt6+3a357dLFTbZKSgqdHskiIiWgJBxsW7a4xHHttdC4MXT3upbjYd55B5o3x06axDtVrmY/MWQbH5ViomjRP87r6Fx97Ndfh4cegjVrQuNDi4hIKehydDDNmOHWrxrj+tYOGJDv8pm5qWkBndRU0OOnL1nG9v116JnzFn+c0JrX707hnD3Jf5l4FVSbNrmJajVrwoQJ7tJ9kJZqiYgEmpJwMFjrEu8pp0CnTvCf/0CjRvlumjcpKm+dbl6hDCibko8HJ11tWcfl368lJjOdd9aeSUrL25k5dTi/7BzBoKGRjBwJMTEe9vrNzHSX6v/1L9ee8emn4ayzvIlFRCRAlIQDKW/Z0d69bh1r69au41EhCiuUURZJeOyCr2m6ZR3TZzxA5exMDLBo9c90nv4AZ50Fb82Hs88u9dOUzuefw223wbp1cOWVcP/9HgckIhIYGhMOlI8+gtNPd43jMzIgJ8evuxVU2rGsSj7+8vteBnz2+sEEnE0Ei/ZfRK0OG1i2LAQS8KRJcN558McfrvzkrFlQr57HQYmIBIaScFnbsQNuvBEuuggqVYLFi+Gll/yueBXoko9Xp62k4+YV5BBBFj4OmCjWXVaTphf/5G1Dpv25HzK6dIFBg+Crr6BHDw8DEhEJPCXhsrZ3r5tl/OCDsHq1m9RUDAEp+ZieDsuXk5UFPxx7D10i3ucC32LGnNqf664dxTdnnBK4kpJF2bbNNVfo2dONnTduDE8+6dZPi4iUcxoTLgvffguvvuomETVq5JYhlTCJlHnJx08+gdtuIzvtZzqdvIUla2rSJr4D2W1X8kpOF+rVjGFMIEtKFiQnx116vv9+d7l+xAh3m5otiEgFoiRcGpmZ7qxt5Eh36fmGG+DUU0t9FpfQMrb0SfGjj2D4cPjiC36v0Zjr9r3B17/U5M03oWfPGIw5v3SPXxpbt8J118HSpa4+9osvupnjIiIVjJJwSS1f7mbwrl7tWuc9+yzEBvlssiALFkC3blhryaYSV+6azEk3d2L6OE9LUh9So4arzPXf/7rxc9V7FpEKSkm4JDIz3ThmdjbMmQP/+IfXETnp6RAdTfpnXxJl8wb8LZNuXsYpL3cq06cqdkGRTz6BZ56B6dNd4Y3Vq4PanlFEJBTpXbA4PvzQjV9GRsLcuW4GbygkYGvhtdegcWOSnlxNr+fiSCf6YK/fU26NK9Onyyv4kbZzP5ZDBUXmpqb9deOdO12xjQsvdFcPtmxxtysBi4goCfvlp5+gVy+3fGbyZHfb2We7y6pe27IFLr4YbryRjVl/o9+90Xwf246tLy/CN3oUZtGiMi85WVhBkYOshdmzoVkz95rddx+sX+/GzEVEBNDl6MLlzeAdNsxd6h092tV+DgUpKTBmDHbhQrKI5IGYZ3l2Tz9GPBbB4MEQGRm4kpN+FRTJznazxY8/3i3ZOuecgMQiIhLOlIQL06+fS8IdO8ILL4TOWVxKCnTqhE1PJ8caejGD3y/4B6tecu11A61ezRjS8knEsdUru8IkvXq5cd/33oMTTnAzx0VE5C/07ni09HQ38apaNTf7uX1718M2FGbwHjgAo0eTs/UHSM8gwlpyiCAxYSOtZwdvmHVI1yZHNJkAaLHzR6Z8MBlWr3AFS+65B+rXD05AIiJhSkn4cIsXu0lEcXHuDLhVK/dVhkrcpnDpUncpfONGPqh1HXE2iigy8EVH0WZoXFBH9w8vKPLbb7sZuuotbvp4BhHVq8GUKe5Di4iIFEkTs8DVe77pJnfZOSfHXU4NgGLNKs6zZw8MGIC94AL+2LaPSyPe56bI/5HyqJt4FZFU9hOv/JHQMpalwzry9b6F3LLoNSKu6gUbNriCJaFw1UBEJAzoTDgpCXr3dktpEhPhoYcgpmyaJRytWG0KU1IgORn278c+/zxTawzkzp2j6fXPqmwYD7Vre9jrd9cu16bxxBNh6FC47DLo1s2bWEREwljFTcLWujO2U05xzeInTHCtBwPI7zaF8+dDz57YrCwyieJmO4WlNfvw1huuOZOn5s2D/v3hzDPdxKvGjd2XiIgUW8W7HJ2RAY895kpNWgsNG7oiHAFOwOBHm0JrYepU6NULe+AAJjsbk53BNef/yLp1HifgvLXSCQlQp46rly0iIqVSsZLwZ5+59aoPPuiqXu3P/8w0UAptU7h1q7uke8MNbK10MulEk4WPiMpRXPrvOKpUCWqoR/r0U1d045133AeYFSugdWsPAxIRKR8qRhLevRvuuMMtN9q1C95+G958E445JqhhJLSMZUzP04mtGYMBYmvGMKbn6SQcb7AtWpCZ/Cn3xzzD3/evYuZtSTBqFL7F3ky8AlzBDXBXCbp0gTVr3Lh5ZKQ38YiIlDPGWhvQJ2jVqpVdsWJFQJ+jSLt3u0RyxRWuilPVqt7Gk2fHDqhTh82b4f1uz/DE/yXQ6PyGvPQS/P3vHsaVmQnjx8Nbb7mzYCVdEZFiMcastNYWuca1/J4Jb9kCAwa4MeDq1V2zhQkTQiMBf/yxq3gVG8vU+1Zx+umQ+NNdJD7fkI8/9jgB511qTkyEBg3gzz89DEZEpHwrf0k4KwvGjYPmzeHVV2HVKnd7GQ2qzk1No/3jSZw07D3aP55U+Brf/EyaBPHxkJREzoEspk/YRseOrrdBv34eNhdKT4fBg6FNG/j1V9ei8c03XflJEREJiPK1RGn5cujb1yXeHj3g2Wfd2VwZySu2kbfWN6/YBhyqIlVoRaz77sNOmACAAXKAJ65ezenTL/G+voXPB4sWuVKdjz+u5CsiEgTl50zYWncq+euvroXe3LllmoCh6BZ+RVXE2vp7NT6oeiX7iSHb+KgUHcUZd8WVeQL2+2z9999h0CBXqCQy0s0ef+EFJWARkSAJ7zNha12yjYuDWrVg5kw47jg3BhwARRXbODpJ1963ixGLJvHx5m4srHE/E199mEaNDDP/lULb9GQXdxnPfPbnbB1rYdYsGDjQJeL4eLf+N0CVwkREJH/hm4S3bnUTr959Fx55BEaMcNWvAqigFn55xTbykvHZP27gppXz6LD5S6IzM1j8bQ9eyIC77zaMGgVVqwau5GSRpTF//BHuvNOt+W3VChYudNWvREQk6MIvCWdlwdNPu6QLbhLWoEFBeer8WvgdLLaBS8bNv1jExHlj8FlLNoY+TGVezctJmefmPJWFwsadiyyNee+98NFHbgnSXXep16+IiIfCb0x46FA3izc+3i07uu++oCWSAott5CbAIV2bcMWGZCJy117nEEGzk77k1bd3lWkCLmzcOb/SmCfv+JGWdrf7ZcIEWLfOJWMlYBERT4Xfu/Ddd8N557nCG8Wc0VTiXr6HSWgZ+9f7rF8P33zDWWcl8LBvNBfSjUpkkl2pEheM6ECHc4v3HIUp6nLz4WfrlbKz6LtsDoOWzuC3jl2Ba6B+/TKLRURESif8knDDhu6rmPyasFRcGRkwZgx29Gh2VW/AGfu6Y31xXHZvEj3rJBMZH0eHMp54VdTl5rx9eefltxn8xr9pun0LaRd1J3bq5DKNQ0RESi/8knAJFauXb1FSUuC112DBAvjuOxbUuZY+O57igksqMXEiNGwYuIlXRU0OA0j47gsSJg6AE06AuXOJvfzygMQiIiKlE35jwiXkdy/foqSkQMeO2BdewH63hSFmLH3MNJ6ZXpd33y3RSXqxFNqJad8+d0OnTm7M96uvQAlYRCRkVZgkXGQvX398+y0kJ2MzMjFAFhG0PjOTDRvgmmuKPURdIvlNDhvXuQEJzz0Mbdu6S+Q1asDYse67iIiErAqThAs9gyzKb79Bnz7YJk2Ysqg++3OiDvb6ver5OI49NkBBFyChZSxLh3Xku8cvZelpO7n06s6uTna3bpCTE9xgRESkxCrMmHDeuG+xZkdbCzNmwKBB5OzcxdPHPEjioqv4vdcp3Nksmcpd47zr9btzJ9xyi2u0cNZZ8N57cPbZ3sQiIiIlUmGSMBSwvKgg1rplUG+9xaY655KQ9TI0bEHyZGjbNnATr/xWpQps2wZjxri10ur5KyISdipUEvbLZ5/Bxx9jO8SRGtmG2cdcyPhdA0kc6SMxEaKiPIxt82YYPhyee87Vyl661MPehyIiUlpKwoebMQOuvx5r4YCpTP+cRZh27fhyMjRr5mFc2dnwzDMuAft8rt1gfLwSsIhImNO7OLgZxaNGYfv0webkYGwOvpwMxndP5pNPPE7A69ZB+/ZuyVF8vKvOFR/vYUAiIlJWlISXLYNzzoERI1hWJZ50osnChy86ivMeiMPnK/ohAmr4cNi0CaZNc52PyrhHsoiIeKfCX47OXLWefd/v5EbfO3wa2Z1pD6fQJSoZEx/n3cznZctcX+TGjWHiRNdooW5db2IREZGAqXhJOCUFXnoJGjXis4se5tYn/8n3u3uRcF1VNjwJdet6OPP5zz9di8annoJrr4WpU+HEE72JRUREAq5iJeEPPoDu3SE7mxwMQ0Z2Zm+D9syaX5Vu3TyOLSnJTbjavBnuuAOeeMLjgEREJNAqxpiwtTBzJlx5JTbbNXHIJoJh7Zawfj3eJ+Bp01y954gISE52l6CrV/c4KBERCbSKkYRTU+Gaa/jR15B0Kh8sOXnZ+DiqVfMwrt273ffLLoNHHoE1a6BDBw8DEhGRYCpXl6PnpqYdLEsZW70yo2P3ceENPZi+4WzmVf+At/d15tmbl/HPxsn4Osd5N/Hq11/hrrvccqOVK91Z74gR3sQiIiKeKTdJeG5qGolz1tJ0yzpu2vgp521dzWm/fc+VT3/JnNQzadOmKysmQ4sWHk68shamT4dBg9xZ8EMPBaf1koiIhKRyk4THLviaFptXMX3mcCrluHHf0RGJzF3XjKefhv79KXLN7+Fn0n41eCiOP/6APn1co4U2beCVVzyuAiIiIl4rN2PCv/y+l2feHkdkTnZur18fB2oYTrzlE+66y78EnDhnLWk792OBtJ37SZyzlrmpaWUTYNWqrvPRhAmu5rMSsIhIhRf+SfjAAQDqVq/GnOqXcoAoMvGRGVGJdZfUolEj/x5m7IKv2Z+ZfcRt+zOzGbvg65LHtmkT9O7tzoIjI2HJErjnnqI/EYiISIUQ3kn4o4+gaVM2jnuXtFfPZ8C2KVzcaBbj297A9dc+ysaTWjCkaxO/Hmrbzv3Fur1Q2dkwfjyccYZbm7x6tbtdDRdEROQw4TcmnJIC8+fDl1/C/Pn8WvM0bhlSGxpE8uDTv7F4X1Ve3HkF9WrGMKYYY7r1asaQlk/CrVczpnjxrV8PN9/sSk927+7W/NavX7zHEBGRCsGvJGyMOR5401p7gTEmEpgD1AZetta+EsgAj5CSExH9xQAACVNJREFUAnFxkJGBBWYfcwN9dr7ILf2j+WAMVKt2LNCxRA89pGsTEuesPeKSdEykz+8z6YOGD3dVr6ZNg2uu0exnEREpUJFJ2BhTC5gCVMm9aSCw0lo70hgz3xgzy1q7J5BBHpScjM3MxADZ+Pi+yt/5aGE07duX/qHzzphLNDt65UqoU8c1XHj+eTVcEBERv/hzJpwN9Abm5f4eBwzL/XkJ0ApYXOaR5ePbBnGcaKOJJANbKYoBs+KIKoMEnCehZWzxliSlp8PIkTBuHFx1lVsDrIYLIiLipyKTsLV2N4A5dFm1CpC3bud34Pij72OM6Qv0BWjYsGFZxAlAg6va8fj0RfQ9NZkTro7zruIVuGVGt9wCX3/txoDHj/cuFhERCUslmZi1F4gBdgFVc38/grV2EjAJoFWrVrY0AR4uKgpGzPew4lWeWbPc0qOGDWHhQrjoIm/jERGRsFSSNTMrgfNzfz4T2FJm0YS6P/9037t2hQcegHXrlIBFRKTESnImPAWYb4y5AGgGfFG2IYWgXbtg8GD4/HNYscI1XHj0Ua+jEhGRMOf3mbC1Ni73+1bgImAp0Nlam13Y/cLeu+9C8+au1nO3bq4Jg4iISBkoUbEOa+024I0yjiW07NkDd94J//ufS8JvvQWtW3sdlYiIlCNhVTEroF2OjhYdDd984/r8PvAAVK4cmOcREZEKK2yScF6Xo7yKVnldjoCyS8S//OJ6/D7+ONSuDZ9+6gpviIiIBEDYdBQISJejPNa6y87NmsFrr7nymKAELCIiARU2SbhMuxwdLi0NevSAPn3gtNMgNRUuvbR0jykiIuKHsEnCBXUzKnaXo6MNHgyLFsGTT7rLz02blu7xRERE/BQ2SXhI1ybERPqOuK1EXY4Atm6F7793P48fD2vWwN13g89X+P1ERETKUNgk4YSWsYzpeTqxNWMwQGzNGMb0PL14k7Jyclx/3xYtYOBAd1u9enDKKQGJWUREpDBhNfOo2F2ODvftt3DrrZCc7EpNPvNMmcYmIiJSXGGVhEssKQm6d4fISJg82XU9OtQVSkRExBNhczm6RLJzlzS1bg3XXgvr17v2g0rAIiISAspnEs7OhrFjoU0bOHAAqlVzZ8D163sdmYiIyEHlLwl/9RWcdx4MHQoNGsC+fV5HJCIikq/yk4SzsuCxx6BlS9i8GWbOhDlzoFYtryMTERHJV/lJwtbC7NmQkODGfnv31tiviIiEtPCeHZ2R4Spd3Xaba7iQnOzGf0VERMJA+J4Jp6a6Wc/DhsGsWe42JWAREQkj4ZeEDxyA4cNdAv71V5g3D26/3euoREREii38kvDQoTB6NFx3nRv77dHD64hERERKJPzGhO+/H7p0UbtBEREJe+GXhOvVc18iIiJhLvwuR4uIiJQTSsIiIiIeURIWERHxiJKwiIiIR5SERUREPKIkLCIi4hElYREREY8oCYuIiHhESVhERMQjSsIiIiIeURIWERHxiJKwiIiIR5SERUREPGKstYF9AmO2A1vL+GGPBX4r48f0QnnZD9C+hKrysi/lZT9A+xKKArEfjay1dYvaKOBJOBCMMSusta28jqO0yst+gPYlVJWXfSkv+wHal1Dk5X7ocrSIiIhHlIRFREQ8Eq5JeJLXAZSR8rIfoH0JVeVlX8rLfoD2JRR5th9hOSYsIiJSHoTrmbCIiEjYC9kkbIx52RiTYowZXpptvGaMqWGMed8Ys9AY85YxJiqfbSoZY743xiTnfp3uRaxF8TdOY8wjxpjlxpjngh2jv4wx/Q7bj1XGmBfz2Sakj4sx5nhjzCe5P0caY94xxiw1xtxcyH382i7YjtqXhrmvd5IxZpIxxhRwn1hjzI+HHZ8il4MEw1H74neMofh+dtS+PHLYfmw0xiQWcJ+QOi75vQf7+1oH45iEZBI2xvQEfNbadsDJxphTS7JNiLgOmGCt7QL8DFyczzZnADOstXG5X2uDGqH/iozTGHMOcD5wLvCrMaZzsIP0h7V2Yt5+AJ8AL+WzWcgeF2NMLWAKUCX3poHASmtte+BKY0y1Au7q73ZBk8++3A70s9Z2BBoABX34aQOMPuz4bA98tIXLZ1/8ijEU38+O3hdr7cOH/c2sA14r4K6hdlyOfg++Gj9e62Adk5BMwkAc8Ebuzwtxb+ol2cZz1trnrbUf5v5aF/g1n83aAt2NMctyP3lVCl6ExeJPnB2A2dZNNlgAXBDUCIvJGBMLHG+tXZHPP4fycckGegO7c3+P49DfwxKgoDWP/m4XTEfsi7X2QWvthtx/q0PBRRTaArcaY740xjwW+DD9cvRx8TfGOELv/ezofQHAGNMa+NFam1bA/ULquOTzHnw9/r3WcX5uVyqhmoSrAHkH+Hfg+BJuEzKMMe2AWtbaz/P55+VAZ2vtuUAkcElQg/OfP3GG1XEB+gMTC/i3kD0u1trd1tpdh93k7+secscnn30BwBjTG1hvrd1WwF3fx71RtgbaGWPOCFyU/slnX/yNMWyOCzAI+E8hdw254wKH3oOBHwihv5VQTcJ7gZjcn6uSf5z+bBMSjDG1cf9pCxqDW2Ot/Sn35xWA55eiCuBPnOF0XCKAeCC5gE3C5biA/697WBwfY8zJwGDg7kI2+8xau8damw2kEprHx98Yw+W41ASOs9Z+W8hmIXdcjnoPDqm/lZA80MBKDp36nwlsKeE2nsudiDULSLTWFlRDe6ox5kxjjA9IAFYHLcDi8SfOsDguuS4AvrAFr9MLl+MC/r/uIX98csciZwA3F3AmlmeBMeZEY8wxQBfcOGWo8TfGkD8uuS4H5hexTUgdl3zeg0Prb8VaG3JfQHXcG94EYEPuC/BoEdvU8DruAvalH/AH7mwrGXg4n31pAawB1uImNHgedwH7ckScQG1g8lHbRABLgaeBr4GTvI67kP15DOiZ+3OzcDwuQHLu90bA+tzXfTngAzoCA47a/i/beb0P+ezLE8BPh/3NdChgX+KBjbnHaECw4/VzX/4SYwH/10L2/SxvX3J/ng6cfdjvIX9c8nkPvvHo19rLYxKyxTpyPw1fBCyx1v5c0m0k+IwxMcClwJfW2s1ex1NRGGPq4T65L7CFnEH6u50El97Pgsff1zoYxyRkk7CIiEh5F6pjwiIiIuWekrCIiIhHlIRFREQ8oiQsIiLiESVhERERjygJi4iIeOT/AS96zYJqzanmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(x, y, 'o', label=\"Data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "ax.plot(x, res2.fittedvalues, 'r--.', label=\"Predicted\")\n",
    "ax.plot(x, iv_u, 'r--')\n",
    "ax.plot(x, iv_l, 'r--')\n",
    "legend = ax.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joint hypothesis test\n",
    "\n",
    "### F test\n",
    "\n",
    "We want to test the hypothesis that both coefficients on the dummy variables are equal to zero, that is, $R \\times \\beta = 0$. An F test leads us to strongly reject the null hypothesis of identical constant in the 3 groups:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 0 0]\n",
      " [0 0 1 0]]\n",
      "<F test: F=array([[145.49268198]]), p=1.2834419617291377e-20, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "R = [[0, 1, 0, 0], [0, 0, 1, 0]]\n",
    "print(np.array(R))\n",
    "print(res2.f_test(R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also use formula-like syntax to test hypotheses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[145.49268198]]), p=1.2834419617291078e-20, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res2.f_test(\"x2 = x3 = 0\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Small group effects\n",
    "\n",
    "If we generate artificial data with smaller group effects, the T test can no longer reject the Null hypothesis: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = [1., 0.3, -0.0, 10]\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + np.random.normal(size=nsample)\n",
    "\n",
    "res3 = sm.OLS(y, X).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.3031864410632063, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res3.f_test(R))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.3031864410632063, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res3.f_test(\"x2 = x3 = 0\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multicollinearity\n",
    "\n",
    "The Longley dataset is well known to have high multicollinearity. That is, the exogenous predictors are highly correlated. This is problematic because it can affect the stability of our coefficient estimates as we make minor changes to model specification. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.datasets.longley import load_pandas\n",
    "y = load_pandas().endog\n",
    "X = load_pandas().exog\n",
    "X = sm.add_constant(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 TOTEMP   R-squared:                       0.995\n",
      "Model:                            OLS   Adj. R-squared:                  0.992\n",
      "Method:                 Least Squares   F-statistic:                     330.3\n",
      "Date:                Tue, 24 Dec 2019   Prob (F-statistic):           4.98e-10\n",
      "Time:                        14:59:01   Log-Likelihood:                -109.62\n",
      "No. Observations:                  16   AIC:                             233.2\n",
      "Df Residuals:                       9   BIC:                             238.6\n",
      "Df Model:                           6                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const      -3.482e+06    8.9e+05     -3.911      0.004    -5.5e+06   -1.47e+06\n",
      "GNPDEFL       15.0619     84.915      0.177      0.863    -177.029     207.153\n",
      "GNP           -0.0358      0.033     -1.070      0.313      -0.112       0.040\n",
      "UNEMP         -2.0202      0.488     -4.136      0.003      -3.125      -0.915\n",
      "ARMED         -1.0332      0.214     -4.822      0.001      -1.518      -0.549\n",
      "POP           -0.0511      0.226     -0.226      0.826      -0.563       0.460\n",
      "YEAR        1829.1515    455.478      4.016      0.003     798.788    2859.515\n",
      "==============================================================================\n",
      "Omnibus:                        0.749   Durbin-Watson:                   2.559\n",
      "Prob(Omnibus):                  0.688   Jarque-Bera (JB):                0.684\n",
      "Skew:                           0.420   Prob(JB):                        0.710\n",
      "Kurtosis:                       2.434   Cond. No.                     4.86e+09\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large, 4.86e+09. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\stats\\stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    }
   ],
   "source": [
    "ols_model = sm.OLS(y, X)\n",
    "ols_results = ols_model.fit()\n",
    "print(ols_results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Condition number\n",
    "\n",
    "One way to assess multicollinearity is to compute the condition number. Values over 20 are worrisome (see Greene 4.9). The first step is to normalize the independent variables to have unit length: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_x = X.values\n",
    "for i, name in enumerate(X):\n",
    "    if name == \"const\":\n",
    "        continue\n",
    "    norm_x[:,i] = X[name]/np.linalg.norm(X[name])\n",
    "norm_xtx = np.dot(norm_x.T,norm_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we take the square root of the ratio of the biggest to the smallest eigen values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56240.8709117813\n"
     ]
    }
   ],
   "source": [
    "eigs = np.linalg.eigvals(norm_xtx)\n",
    "condition_number = np.sqrt(eigs.max() / eigs.min())\n",
    "print(condition_number)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dropping an observation\n",
    "\n",
    "Greene also points out that dropping a single observation can have a dramatic effect on the coefficient estimates: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage change 4.55%\n",
      "Percentage change -2228.01%\n",
      "Percentage change 154304695.31%\n",
      "Percentage change 1366329.02%\n",
      "Percentage change 1112549.36%\n",
      "Percentage change 92708715.91%\n",
      "Percentage change 817944.26%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ols_results2 = sm.OLS(y.iloc[:14], X.iloc[:14]).fit()\n",
    "print(\"Percentage change %4.2f%%\\n\"*7 % tuple([i for i in (ols_results2.params - ols_results.params)/ols_results.params*100]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at formal statistics for this such as the DFBETAS -- a standardized measure of how much each coefficient changes when that observation is left out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "infl = ols_results.get_influence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general we may consider DBETAS in absolute value greater than $2/\\sqrt{N}$ to be influential observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2./len(X)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\stats\\outliers_influence.py:695: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return self.resid / sigma / np.sqrt(1 - hii)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    dfb_const  dfb_GNPDEFL       dfb_GNP     dfb_UNEMP     dfb_ARMED  \\\n",
      "0   -0.016406  -169.822675  1.673981e+06  54490.318088  51447.824036   \n",
      "1   -0.020608  -187.251727  1.829990e+06  54495.312977  52659.808664   \n",
      "2   -0.008382   -65.417834  1.587601e+06  52002.330476  49078.352378   \n",
      "3    0.018093   288.503914  1.155359e+06  56211.331922  60350.723082   \n",
      "4    1.871260  -171.109595  4.498197e+06  82532.785818  71034.429294   \n",
      "5   -0.321373  -104.123822  1.398891e+06  52559.760056  47486.527649   \n",
      "6    0.315945  -169.413317  2.364827e+06  59754.651394  50371.817827   \n",
      "7    0.015816   -69.343793  1.641243e+06  51849.056936  48628.749338   \n",
      "8   -0.004019   -86.903523  1.649443e+06  52023.265116  49114.178265   \n",
      "9   -1.018242  -201.315802  1.371257e+06  56432.027292  53997.742487   \n",
      "10   0.030947   -78.359439  1.658753e+06  52254.848135  49341.055289   \n",
      "11   0.005987  -100.926843  1.662425e+06  51744.606934  48968.560299   \n",
      "12  -0.135883   -32.093127  1.245487e+06  50203.467593  51148.376274   \n",
      "13   0.032736   -78.513866  1.648417e+06  52509.194459  50212.844641   \n",
      "14   0.305868   -16.833121  1.829996e+06  60975.868083  58263.878679   \n",
      "15  -0.538323   102.027105  1.344844e+06  54721.897640  49660.474568   \n",
      "\n",
      "          dfb_POP      dfb_YEAR  \n",
      "0   207954.113590 -31969.158503  \n",
      "1    25343.938291 -29760.155888  \n",
      "2   107465.770565 -29593.195253  \n",
      "3   456190.215133 -36213.129569  \n",
      "4  -389122.401699 -49905.782854  \n",
      "5   144354.586054 -28985.057609  \n",
      "6  -107413.074918 -32984.462465  \n",
      "7    92843.959345 -29724.975873  \n",
      "8    83931.635336 -29563.619222  \n",
      "9    18392.575057 -29203.217108  \n",
      "10   93617.648517 -29846.022426  \n",
      "11   95414.217290 -29690.904188  \n",
      "12  258559.048569 -29296.334617  \n",
      "13  104434.061226 -30025.564763  \n",
      "14  275103.677859 -36060.612522  \n",
      "15 -110176.960671 -28053.834556  \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\stats\\_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in greater\n",
      "  return (self.a < x) & (x < self.b)\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\stats\\_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in less\n",
      "  return (self.a < x) & (x < self.b)\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\stats\\_distn_infrastructure.py:1821: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= self.a)\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\stats\\outliers_influence.py:729: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_internal * np.sqrt(hii / (1 - hii))\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\stats\\outliers_influence.py:758: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_external * np.sqrt(hii / (1 - hii))\n"
     ]
    }
   ],
   "source": [
    "print(infl.summary_frame().filter(regex=\"dfb\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
